CoastalME (Coastal Modelling Environment)
Simulates the long-term behaviour of complex coastlines
Loading...
Searching...
No Matches
do_cliff_collapse.cpp
Go to the documentation of this file.
1
10
11/* ==============================================================================================================================
12 This file is part of CoastalME, the Coastal Modelling Environment.
13
14 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.
15
16 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.
17
18 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.
19==============================================================================================================================*/
20#include <cstddef>
21#include <assert.h>
22
23#include <iostream>
24using std::cerr;
25using std::endl;
26using std::ios;
27
28#include "cme.h"
29#include "simulation.h"
30#include "cliff.h"
31#include "cell_talus.h"
32#include "coast_landform.h"
33#include "2di_point.h"
34
35//===============================================================================================================================
38//===============================================================================================================================
40{
42 LogStream << m_ulIter << ": Calculating cliff collapse" << endl;
43
44 int nRet;
45
46 // First go along each coastline and at each point on the coastline, update the total wave energy which it has experienced TODO Note that currently, only cliff objects respond to accumulated wave energy
47 for (int nCoast = 0; nCoast < static_cast<int>(m_VCoast.size()); nCoast++)
48 {
49 for (int nCoastPoint = 0; nCoastPoint < m_VCoast[nCoast].nGetCoastlineSize(); nCoastPoint++)
50 {
51 CACoastLandform* pCoastLandform = m_VCoast[nCoast].pGetCoastLandform(nCoastPoint);
52
53 // Get the coords of the grid cell marked as coastline for the coastal landform object
54 int const nX = m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nCoastPoint)->nGetX();
55 int const nY = m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nCoastPoint)->nGetY();
56
57 // TODO Talus is protecting this cell
58 if (m_pRasterGrid->m_Cell[nX][nY].dGetTalusDepth() > 0)
59 continue;
60
61 // First get wave energy for the coastal landform object
62 double const dWaveHeightAtCoast = m_VCoast[nCoast].dGetCoastWaveHeight(nCoastPoint);
63
64 // If the waves at this point are off-shore, then do nothing, just move to next coast point
65 if (bFPIsEqual(dWaveHeightAtCoast, DBL_NODATA, TOLERANCE))
66 continue;
67
68 // OK we have on-shore waves so get the previously-calculated wave energy
69 double const dWaveEnergy = m_VCoast[nCoast].dGetWaveEnergyAtBreaking(nCoastPoint);
70
71 // And save the accumulated value
72 pCoastLandform->IncTotAccumWaveEnergy(dWaveEnergy);
73
74 int const nCat = pCoastLandform->nGetLandFormCategory();
75
76 // Is this a cliff?
77 if ((nCat == LF_CLIFF_ON_COASTLINE) || (nCat == LF_CLIFF_INLAND))
78 {
79 // It is, so get the cliff object
80 CRWCliff* pCliff = reinterpret_cast<CRWCliff*>(pCoastLandform);
81
82 // And do the notch incision, if any. Note that we consider sediment eroded due to notch incision to be still in place until cliff collapse, i.e. the sediment which filled the notch, pre-incision, is assumed to remain there. If the notch is eventually incised sufficiently to cause cliff collapse, then the sediment from the notch volume is included with the above-notch talus
83 if (! bIncreaseCliffNotchIncision(nCoast, nX, nY, pCliff, dWaveEnergy))
84 // No incision of this cliff
85 continue;
86
87 // OK, we've had some incision. So is the notch now extended enough to cause collapse (either because the overhang is greater than the threshold overhang, or because there is no sediment remaining)?
89 {
90 // // DEBUG CODE ============================================================================================================================================
91 // // Get total depths of sand consolidated and unconsolidated for every cell
92 // if (m_ulIter == 5)
93 // {
94 // double dTmpSandCons = 0;
95 // double dTmpSandUncons = 0;
96 // for (int nX1 = 0; nX1 < m_nXGridSize; nX1++)
97 // {
98 // for (int nY1 = 0; nY1 < m_nYGridSize; nY1++)
99 // {
100 // dTmpSandCons += m_pRasterGrid->m_Cell[nX1][nY1].dGetConsSandDepthAllLayers();
101 //
102 // dTmpSandUncons += m_pRasterGrid->m_Cell[nX1][nY1].dGetUnconsSandDepthAllLayers();
103 // }
104 // }
105 //
106 // // Get the cliff cell's grid coords
107 // int nXCliff = pCliff->pPtiGetCellMarkedAsCliff()->nGetX();
108 // int nYCliff = pCliff->pPtiGetCellMarkedAsCliff()->nGetY();
109 //
110 // // Get this cell's polygon
111 // int nPoly = m_pRasterGrid->m_Cell[nXCliff][nYCliff].nGetPolygonID();
112 //
113 // LogStream << endl;
114 // LogStream << "*****************************" << endl;
115 // LogStream << m_ulIter << ": before cliff collapse on nPoly = " << nPoly << " total consolidated sand = " << dTmpSandCons * m_dCellArea << " total unconsolidated sand = " << dTmpSandUncons * m_dCellArea << endl;
116 // }
117 // // DEBUG CODE ============================================================================================================================================
118
119 // It is ready to collapse
120 int nNotchLayer;
121 double dCliffElevPreCollapse = 0;
122 double dCliffElevPostCollapse = 0;
123 double dFineCollapse = 0;
124 double dSandCollapse = 0;
125 double dCoarseCollapse = 0;
126
127 // Do the cliff collapse
128 nRet = nDoCliffCollapse(nCoast, pCliff, dFineCollapse, dSandCollapse, dCoarseCollapse, nNotchLayer, dCliffElevPreCollapse, dCliffElevPostCollapse);
129 if (nRet != RTN_OK)
130 {
132 LogStream << m_ulIter << ": " << WARN << "problem with cliff collapse, continuing however" << endl;
133 }
134
135 // Deposit all sand and/or coarse sediment derived from this cliff collapse as talus, on the cell on which collapse occurred
136 DoCliffCollapseTalusDeposition(nCoast, pCliff, dSandCollapse, dCoarseCollapse, nNotchLayer);
137
138 // // DEBUG CODE ============================================================================================================================================
139 // // Get total depths of sand consolidated and unconsolidated for every cell
140 // if (m_ulIter == 5)
141 // {
142 // double dTmpSandCons = 0;
143 // double dTmpSandUncons = 0;
144 // for (int nX1 = 0; nX1 < m_nXGridSize; nX1++)
145 // {
146 // for (int nY1 = 0; nY1 < m_nYGridSize; nY1++)
147 // {
148 // dTmpSandCons += m_pRasterGrid->m_Cell[nX1][nY1].dGetConsSandDepthAllLayers();
149 //
150 // dTmpSandUncons += m_pRasterGrid->m_Cell[nX1][nY1].dGetUnconsSandDepthAllLayers();
151 // }
152 // }
153 //
154 // // Get the cliff cell's grid coords
155 // int nXCliff = pCliff->pPtiGetCellMarkedAsCliff()->nGetX();
156 // int nYCliff = pCliff->pPtiGetCellMarkedAsCliff()->nGetY();
157 //
158 // // Get this cell's polygon
159 // int nPoly = m_pRasterGrid->m_Cell[nXCliff][nYCliff].nGetPolygonID();
160 //
161 // LogStream << endl;
162 // LogStream << "*****************************" << endl;
163 // LogStream << m_ulIter << ": after cliff collapse on nPoly = " << nPoly << " total consolidated sand = " << dTmpSandCons * m_dCellArea << " total unconsolidated sand = " << dTmpSandUncons * m_dCellArea << endl;
164 // LogStream << m_ulIter << ": total consolidated sand lost this iteration = " << (m_dStartIterConsSandAllCells - dTmpSandCons) * m_dCellArea << endl;
165 // LogStream << m_ulIter << ": total unconsolidated sand added this iteration = " << (dTmpSandUncons - m_dStartIterUnconsSandAllCells) * m_dCellArea << endl;
166 //
167 // double dTmpAllPolySandErosion = 0;
168 // double dTmpAllPolySandDeposition = 0;
169 // for (unsigned int n = 0; n < m_pVCoastPolygon.size(); n++)
170 // {
171 // double dTmpSandErosion = m_pVCoastPolygon[n]->dGetCliffCollapseErosionSand() * m_dCellArea ;
172 // double dTmpSandDeposition = m_pVCoastPolygon[n]->dGetCliffCollapseUnconsSandDeposition() * m_dCellArea ;
173 //
174 // LogStream << m_ulIter << ": polygon = " << m_pVCoastPolygon[n]->nGetPolygonCoastID() << " sand erosion = " << dTmpSandErosion << " sand deposition = " << dTmpSandDeposition << endl;
175 //
176 // dTmpAllPolySandErosion += dTmpSandErosion;
177 // dTmpAllPolySandDeposition += dTmpSandDeposition;
178 // }
179 //
180 // LogStream << "-------------------------------------------" << endl;
181 // LogStream << m_ulIter << ": all polygons, sand erosion = " << dTmpAllPolySandErosion << " sand deposition = " << dTmpAllPolySandDeposition << endl;
182 // LogStream << "*****************************" << endl;
183 // }
184 // // DEBUG CODE ============================================================================================================================================
185 }
186 }
187 }
188 }
189
192
193 return RTN_OK;
194}
195
196//===============================================================================================================================
198//===============================================================================================================================
199int CSimulation::nDoCliffCollapse(int const nCoast, CRWCliff* pCliff, double& dFineCollapse, double& dSandCollapse, double& dCoarseCollapse, int& nNotchLayer, double& dPreCollapseCellElev, double& dPostCollapseCellElevNoTalus)
200{
201 // Get the cliff cell's grid coords
202 int const nX = pCliff->pPtiGetCellMarkedAsCliff()->nGetX();
203 int const nY = pCliff->pPtiGetCellMarkedAsCliff()->nGetY();
204
205 // Get this cell's polygon
206 int const nPoly = m_pRasterGrid->m_Cell[nX][nY].nGetPolygonID();
207 if (nPoly == INT_NODATA)
208 {
209 // This cell isn't in a polygon
210 LogStream << m_ulIter << ": " << WARN << "in nDoCliffCollapse(), [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "} is not in a polygon" << endl;
212 }
213
214 // Get a pointer to the polygon
215 CGeomCoastPolygon* pPolygon = m_VCoast[nCoast].pGetPolygon(nPoly);
216
217 // Get the elevation of the apex of the notch from the cliff object
218 double const dNotchElev = pCliff->dGetNotchApexElev();
219
220 // Get the index of the layer containing the notch (layer 0 being just above basement)
221 nNotchLayer = m_pRasterGrid->m_Cell[nX][nY].nGetLayerAtElev(dNotchElev);
222
223 // Safety check: is the notch elevation above the top of the consolidated sediment? If so, do no more
224 if (nNotchLayer == ELEV_ABOVE_SEDIMENT_TOP)
225 {
226#ifdef _DEBUG
227 double const dTopElevNoTalus = m_pRasterGrid->m_Cell[nX][nY].dGetConsSedTopElevOmitTalus();
228 double const dTopElevIncTalus = m_pRasterGrid->m_Cell[nX][nY].dGetConsSedTopElevIncTalus();
229 LogStream << m_ulIter << ": cliff ready to collapse at [" << nX << "][" << nY << "] but notch apex is above sediment top, nNotchLayer = " << nNotchLayer << " dNotchElev = " << dNotchElev << " cons sediment top elev without talus = " << dTopElevNoTalus << " cons sediment top elev inc talus = " << dTopElevIncTalus << endl;
230#endif
231 return RTN_OK;
232 }
233
234 // More safety checks
235 if (nNotchLayer == ELEV_IN_BASEMENT)
236 {
237 LogStream << m_ulIter << ": " << WARN << "in nDoCliffCollapse(), [" << nX << "][" << nY << "] nNotchLayer is in basement" << endl;
238 return RTN_ERR_CLIFF_NOTCH;
239 }
240
241 if (nNotchLayer < 0)
242 {
243 LogStream << m_ulIter << ": " << WARN << "in nDoCliffCollapse(), [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "} notch layer = " << nNotchLayer << ", dNotchElev = " << dNotchElev << " m_dNotchApexAboveMHW = " << m_dNotchApexAboveMHW << " dPreCollapseCellElev = " << dPreCollapseCellElev << endl;
244 return RTN_ERR_CLIFF_NOTCH;
245 }
246
247 // Notch layer is OK, so flag the coastline cliff object as having collapsed
248 pCliff->SetCliffCollapsed();
249
250 int const nTopLayer = m_pRasterGrid->m_Cell[nX][nY].nGetNumOfTopLayerAboveBasement();
251
252 // Safety check
253 if (nTopLayer == INT_NODATA)
254 {
255 LogStream << m_ulIter << ": " << WARN << "in nDoCliffCollapse(), [" << nX << "][" << nY << "] nTopLayer = " << nTopLayer << endl;
257 }
258
259 // Set flags to say that the notch layer, and all layers above it, have changed
260 for (int nLayer = nNotchLayer; nLayer <= nTopLayer; nLayer++)
261 {
262 m_bConsChangedThisIter[nLayer] = true;
263 m_bUnconsChangedThisIter[nLayer] = true;
264 }
265
266 // Get the pre-collapse cliff elevation (we assume that there is no talus on top of the cliff)
267 dPreCollapseCellElev = m_pRasterGrid->m_Cell[nX][nY].dGetAllSedTopElevOmitTalus();
268
269 // Now calculate the vertical depth of sediment lost in this cliff collapse, note that this includes the sediment which filled the notch before any incision took place
270 double dAvailable = 0;
271 double dFineConsLost = 0;
272 double dFineUnconsLost = 0;
273 double dSandConsLost = 0;
274 double dSandUnconsLost = 0;
275 double dCoarseConsLost = 0;
276 double dCoarseUnconsLost = 0;
277
278 // Update the cell's sediment. If there are sediment layers above the notched layer, we must remove sediment from the whole depth of each layer
279 for (int n = nTopLayer; n > nNotchLayer; n--)
280 {
281 // Start with the unconsolidated sediment
282 dAvailable = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetUnconsolidatedSediment()->dGetFineDepth() - m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetUnconsolidatedSediment()->dGetNotchFineLost();
283
284 if (dAvailable > 0)
285 {
286 dFineCollapse += dAvailable;
287 dFineUnconsLost += dAvailable;
288 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetUnconsolidatedSediment()->SetFineDepth(0);
289 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetUnconsolidatedSediment()->SetNotchFineLost(0);
290 }
291
292 dAvailable = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetUnconsolidatedSediment()->dGetSandDepth() - m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetUnconsolidatedSediment()->dGetNotchSandLost();
293
294 if (dAvailable > 0)
295 {
296 dSandCollapse += dAvailable;
297 dSandUnconsLost += dAvailable;
298 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetUnconsolidatedSediment()->SetSandDepth(0);
299 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetUnconsolidatedSediment()->SetNotchSandLost(0);
300 }
301
302 dAvailable = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetUnconsolidatedSediment()->dGetCoarseDepth() - m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetUnconsolidatedSediment()->dGetNotchCoarseLost();
303
304 if (dAvailable > 0)
305 {
306 dCoarseCollapse += dAvailable;
307 dCoarseUnconsLost += dAvailable;
308 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetUnconsolidatedSediment()->SetCoarseDepth(0);
309 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetUnconsolidatedSediment()->SetNotchCoarseLost(0);
310 }
311
312 // Now get the consolidated sediment
313 dAvailable = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetConsolidatedSediment()->dGetFineDepth() - m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetConsolidatedSediment()->dGetNotchFineLost();
314
315 if (dAvailable > 0)
316 {
317 dFineCollapse += dAvailable;
318 dFineConsLost += dAvailable;
319 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetConsolidatedSediment()->SetFineDepth(0);
320 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetConsolidatedSediment()->SetNotchFineLost(0);
321 }
322
323 dAvailable = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetConsolidatedSediment()->dGetSandDepth() - m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetConsolidatedSediment()->dGetNotchSandLost();
324
325 if (dAvailable > 0)
326 {
327 dSandCollapse += dAvailable;
328 dSandConsLost += dAvailable;
329 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetConsolidatedSediment()->SetSandDepth(0);
330 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetConsolidatedSediment()->SetNotchSandLost(0);
331 }
332
333 dAvailable = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetConsolidatedSediment()->dGetCoarseDepth() - m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetConsolidatedSediment()->dGetNotchCoarseLost();
334
335 if (dAvailable > 0)
336 {
337 dCoarseCollapse += dAvailable;
338 dCoarseConsLost += dAvailable;
339 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetConsolidatedSediment()->SetCoarseDepth(0);
340 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(n)->pGetConsolidatedSediment()->SetNotchCoarseLost(0);
341 }
342 }
343
344 // Now calculate the sediment lost from the consolidated layer into which the erosional notch was incised
345 double const dNotchLayerTop = m_pRasterGrid->m_Cell[nX][nY].dCalcLayerElev(nNotchLayer);
346 double const dNotchLayerThickness = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->dGetTotalThickness();
347 double const dNotchLayerFracRemoved = (dNotchLayerTop - dNotchElev) / dNotchLayerThickness;
348
349 // Sort out the notched layer's sediment, both consolidated and unconsolidated, for this cell. First the unconsolidated sediment
350 double dFineDepth = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetUnconsolidatedSediment()->dGetFineDepth();
351 dAvailable = dFineDepth - m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetUnconsolidatedSediment()->dGetNotchFineLost();
352
353 if (dAvailable > 0)
354 {
355 // Some unconsolidated fine sediment is available for collapse
356 double const dLost = dAvailable * dNotchLayerFracRemoved;
357 dFineCollapse += dLost;
358 dFineUnconsLost += dLost;
359 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetUnconsolidatedSediment()->SetFineDepth(dFineDepth - dLost);
360 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetUnconsolidatedSediment()->SetNotchFineLost(0);
361 }
362
363 double dSandDepth = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetUnconsolidatedSediment()->dGetSandDepth();
364 dAvailable = dSandDepth - m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetUnconsolidatedSediment()->dGetNotchSandLost();
365
366 if (dAvailable > 0)
367 {
368 // Some unconsolidated sand sediment is available for collapse
369 double const dLost = dAvailable * dNotchLayerFracRemoved;
370 dSandCollapse += dLost;
371 dSandUnconsLost += dLost;
372 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetUnconsolidatedSediment()->SetSandDepth(dSandDepth - dLost);
373 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetUnconsolidatedSediment()->SetNotchSandLost(0);
374 }
375
376 double dCoarseDepth = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetUnconsolidatedSediment()->dGetCoarseDepth();
377 dAvailable = dCoarseDepth - m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetUnconsolidatedSediment()->dGetNotchCoarseLost();
378
379 if (dAvailable > 0)
380 {
381 // Some unconsolidated coarse sediment is available for collapse
382 double const dLost = dAvailable * dNotchLayerFracRemoved;
383 dCoarseCollapse += dLost;
384 dCoarseUnconsLost += dLost;
385 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetUnconsolidatedSediment()->SetCoarseDepth(dCoarseDepth - dLost);
386 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetUnconsolidatedSediment()->SetNotchCoarseLost(0);
387 }
388
389 // Do the same for the consolidated sediment
390 dFineDepth = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetConsolidatedSediment()->dGetFineDepth();
391 dAvailable = dFineDepth - m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetConsolidatedSediment()->dGetNotchFineLost();
392
393 if (dAvailable > 0)
394 {
395 // Some consolidated fine sediment is available for collapse
396 double const dLost = dAvailable * dNotchLayerFracRemoved;
397 dFineCollapse += dLost;
398 dFineConsLost += dLost;
399 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetConsolidatedSediment()->SetFineDepth(dFineDepth - dLost);
400 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetConsolidatedSediment()->SetNotchFineLost(0);
401 }
402
403 dSandDepth = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetConsolidatedSediment()->dGetSandDepth();
404 dAvailable = dSandDepth - m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetConsolidatedSediment()->dGetNotchSandLost();
405
406 if (dAvailable > 0)
407 {
408 // Some consolidated sand sediment is available for collapse
409 double const dLost = dAvailable * dNotchLayerFracRemoved;
410 dSandCollapse += dLost;
411 dSandConsLost += dLost;
412 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetConsolidatedSediment()->SetSandDepth(dSandDepth - dLost);
413 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetConsolidatedSediment()->SetNotchSandLost(0);
414 }
415
416 dCoarseDepth = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetConsolidatedSediment()->dGetCoarseDepth();
417 dAvailable = dCoarseDepth - m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetConsolidatedSediment()->dGetNotchCoarseLost();
418
419 if (dAvailable > 0)
420 {
421 // Some consolidated coarse sediment is available for collapse
422 double const dLost = dAvailable * dNotchLayerFracRemoved;
423 dCoarseCollapse += dLost;
424 dCoarseConsLost += dLost;
425 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetConsolidatedSediment()->SetCoarseDepth(dCoarseDepth - dLost);
426 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer)->pGetConsolidatedSediment()->SetNotchCoarseLost(0);
427 }
428
429 // Update the cell's totals for cliff collapse erosion
430 m_pRasterGrid->m_Cell[nX][nY].IncrCliffCollapseErosion(dFineCollapse, dSandCollapse, dCoarseCollapse);
431
432 // Update the cell's layer elevations (pre talus deposition) and d50
433 m_pRasterGrid->m_Cell[nX][nY].CalcAllLayerElevsAndD50();
434
435 // Get the post-collapse cell elevation (talus will be deposited above this)
436 dPostCollapseCellElevNoTalus = m_pRasterGrid->m_Cell[nX][nY].dGetAllSedTopElevOmitTalus();
437
439 LogStream << m_ulIter << ": \tcoast " << nCoast << " cliff collapse at [" << nX << "][" << nY << "], before talus deposition: original cell elevation = " << dPreCollapseCellElev << ", new cell elevation = " << dPostCollapseCellElevNoTalus << ", change in elevation = " << dPreCollapseCellElev - dPostCollapseCellElevNoTalus << endl;
440
441 // Update this-polygon totals: add to the depths of cliff collapse erosion for this polygon
442 pPolygon->AddCliffCollapseErosionFine(dFineCollapse);
443 pPolygon->AddCliffCollapseToSuspensionFine(dFineCollapse);
444 pPolygon->AddCliffCollapseErosionSand(dSandCollapse);
445 pPolygon->AddCliffCollapseErosionCoarse(dCoarseCollapse);
446
447 // And update the this-timestep totals and the grand totals for the number of cells with cliff collapse
450
451 // Add to this-iteration totals of fine sediment (consolidated and unconsolidated) eroded via cliff collapse
454
455 // Also add to the total suspended load. Note that this addition to the suspended load has not yet been added to all cells, this happens in nUpdateGrid()
456 m_dThisIterFineSedimentToSuspension += (dFineConsLost + dFineUnconsLost);
457
458 // Add to this-iteration totals of sand and coarse sediment (consolidated and unconsolidated) eroded via cliff collapse
463
464#ifdef _DEBUG
465 // Save the timestep at which cliff collapse occurred
466 m_pRasterGrid->m_Cell[nX][nY].pGetLandform()->SetCliffCollapseTimestep(m_ulIter);
467#endif
468
469 // Reset cell cliff info
470 m_pRasterGrid->m_Cell[nX][nY].pGetLandform()->SetCliffNotchIncisionDepth(m_dCellSide);
471
472 // Final safety check
473 int const nNewTopLayer = m_pRasterGrid->m_Cell[nX][nY].nGetNumOfTopLayerAboveBasement();
474 if (nNewTopLayer == INT_NODATA)
476
477 return RTN_OK;
478}
479
480//===============================================================================================================================
482//===============================================================================================================================
483bool CSimulation::bIncreaseCliffNotchIncision(int const nCoast, int const nX, int const nY, CRWCliff* pCliff, double const dWaveEnergy)
484{
485 // Get the coastline point of the cliff
486 int const nCoastPoint = pCliff->nGetPointOnCoast();
487
488 // And get the wave runup at this point
489 double const dRunup = m_VCoast[nCoast].dGetRunUp(nCoastPoint);
490 double const dWaveElev = m_dThisIterSWL + dRunup;
491
492 // Get the apex elevation of the cliff notch
493 double const dNotchApexElev = pCliff->dGetNotchApexElev();
494
495 // Is there a notch?
496 if (! bFPIsEqual(dNotchApexElev, DBL_NODATA, TOLERANCE))
497 {
498 // This is a notch in this cliff object
499 double const dSedTopElevNoTalus = m_pRasterGrid->m_Cell[nX][nY].dGetAllSedTopElevOmitTalus();
500
501 assert(dNotchApexElev <= dSedTopElevNoTalus);
502
503 // Get the cutoff elevation (if this-iteration SWL is below this, there is no incision)
504 double const dCutoffElev = dNotchApexElev - CLIFF_NOTCH_CUTOFF_DISTANCE;
505
506 if (dWaveElev < dCutoffElev)
507 {
508 // SWL is below the cutoff elevation, so no incision of this existing notch
509#ifdef _DEBUG
510 double const dSedTopElevIncTalus = m_pRasterGrid->m_Cell[nX][nY].dGetAllSedTopElevIncTalus();
511
512 LogStream << m_ulIter << ": \tNO incision of existing notch at [" << nX << "][" << nY << "] dWaveElev = " << dWaveElev << " dCutoffElev = " << dCutoffElev << " dRunup = " << dRunup << " dNotchApexElev = " << dNotchApexElev << " dSedTopElevNoTalus = " << dSedTopElevNoTalus << " dSedTopElevIncTalus = " << dSedTopElevIncTalus << endl;
513#endif
514 return false;
515 }
516
517 // SWL is above the cutoff, so we have some incision of this existing notch
518 double dWeight;
519 if (dWaveElev > dNotchApexElev)
520 // Should not happen: safety check
521 dWeight = 1;
522 else
523 // Assume a linear decrease in incision with distance downwards from notch apex
524 dWeight = 1 - ((dNotchApexElev - dWaveElev) / CLIFF_NOTCH_CUTOFF_DISTANCE);
525
526 // Calculate this-timestep cliff notch erosion (is a length in external CRS units). Note that only consolidated sediment can have a cliff notch
527 double const dNotchIncision = dWeight * dWaveEnergy / m_dCliffErosionResistance;
528
529 // Deepen the cliff object's erosional notch as a result of wave energy during this timestep. Note that notch deepening may be constrained, since this-timestep notch extension cannot exceed the length (i.e. cellside minus notch depth) of sediment remaining on the cell
530 pCliff->IncreaseNotchIncision(dNotchIncision);
531
532 // And add to the cell's accumulated wave energy
533 m_pRasterGrid->m_Cell[nX][nY].pGetLandform()->AddToAccumWaveEnergy(dWaveEnergy* dWeight);
534
535#ifdef _DEBUG
536 LogStream << m_ulIter << ": \tincision of existing notch at [" << nX << "][" << nY << "] dWaveElev = " << dWaveElev << " dCutoffElev = " << dCutoffElev << " dRunup = " << dRunup << " dWeight = " << dWeight << " dNotchApexElev = " << dNotchApexElev << " dSedTopElevNoTalus = " << dSedTopElevNoTalus << " dNotchIncision = " << dNotchIncision << endl;
537#endif
538 return true;
539 }
540 else
541 {
542 // No notch in this cliff object. Can we create one?
543 double const dSedTopElevNoTalus = m_pRasterGrid->m_Cell[nX][nY].dGetAllSedTopElevOmitTalus();
544
545 if (m_dThisIterNewNotchApexElev < dSedTopElevNoTalus)
546 {
547 // Yes we can create a notch here
549 pCliff->SetNotchIncision(0);
550
551 // Get the cutoff elevation (if this-iteration SWL is below this, there is no incision)
552 double const dCutoffElev = m_dThisIterNewNotchApexElev - CLIFF_NOTCH_CUTOFF_DISTANCE;
553
554 if (dWaveElev < dCutoffElev)
555 {
556 // SWL is below the cutoff elevation, so no incision of this existing notch
557#ifdef _DEBUG
558 double const dSedTopElevIncTalus = m_pRasterGrid->m_Cell[nX][nY].dGetAllSedTopElevIncTalus();
559
560 LogStream << m_ulIter << ": \tNO incision of new notch at [" << nX << "][" << nY << "] dWaveElev = " << dWaveElev << " dCutoffElev = " << dCutoffElev << " dRunup = " << dRunup << " dNotchApexElev = " << dNotchApexElev << " dSedTopElevNoTalus = " << dSedTopElevNoTalus << " dSedTopElevIncTalus = " << dSedTopElevIncTalus << endl;
561#endif
562 return false;
563 }
564
565 // We have some notch incision of this newly-created notch
566 double dWeight;
567 if (dWaveElev > dNotchApexElev)
568 // Should not happen: safety check
569 dWeight = 1;
570 else
571 // Assume a linear decrease in incision with distance downwards from notch apex
572 dWeight = 1 - ((dNotchApexElev - dWaveElev) / CLIFF_NOTCH_CUTOFF_DISTANCE);
573
574 // Calculate this-timestep cliff notch erosion (is a length in external CRS units). Note that only consolidated sediment can have a cliff notch
575 double const dNotchIncision = dWeight * dWaveEnergy / m_dCliffErosionResistance;
576
577 // Deepen the cliff object's erosional notch as a result of wave energy during this timestep. Note that notch deepening may be constrained, since this-timestep notch extension cannot exceed the length (i.e. cellside minus notch depth) of sediment remaining on the cell
578 pCliff->IncreaseNotchIncision(dNotchIncision);
579
580 // And add to the cell's accumulated wave energy
581 m_pRasterGrid->m_Cell[nX][nY].pGetLandform()->AddToAccumWaveEnergy(dWaveEnergy * dWeight);
582
583#ifdef _DEBUG
584 LogStream << m_ulIter << ": \tincision of newly-created notch at [" << nX << "][" << nY << "] dWaveElev = " << dWaveElev << " dCutoffElev = " << dCutoffElev << " dRunup = " << dRunup << " dWeight = " << dWeight << " dNotchApexElev = " << dNotchApexElev << " dSedTopElevNoTalus = " << dSedTopElevNoTalus << " dNotchIncision = " << dNotchIncision << endl;
585#endif
586 return true;
587 }
588 else
589 {
590 // The top of the notch would be above the top of the sediment on this cell, so we must try to create a notch further inland. Can't do this for points at the beginning and the end of the coast however
591 if ((nCoastPoint == 0) || (nCoastPoint == m_VCoast[nCoast].nGetCoastlineSize()-1))
592 return false;
593
594 return bCreateNotchInland(nCoast, nCoastPoint, nX, nY, dWaveEnergy, dWaveElev);
595 }
596 }
597
598 // Should never get here
599 return false;
600}
601
602//===============================================================================================================================
604//===============================================================================================================================
605bool CSimulation::bCreateNotchInland(int const nCoast, int const nCoastPoint, int const nX, int const nY, double const dWaveEnergy, double const dWaveElev)
606{
607 LogStream << m_ulIter << ": \tIn bCreateNotchInland() for [" << nX << "][" << nY << "]" << endl;
608
609 bool bFound = false;
610 int const nSeaHandedness = m_VCoast[nCoast].nGetSeaHandedness();
611
612 // Get the points before and after this coastline point
613 CGeom2DIPoint const* pPtiBefore = m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nCoastPoint-1);
614 CGeom2DIPoint const* pPtiAfter = m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nCoastPoint+1);
615
616 // Loop until we can create an inland notch, or until we hit the edge of the grid
617 int n = 1;
618 do
619 {
620 // Get an inland point, planview orthogonal to the coastline at the coast point
621 CGeom2DIPoint const PtiTmp = PtiGetPerpendicular(pPtiBefore, pPtiAfter, n, nSeaHandedness);
622
623 if (! bIsWithinValidGrid(&PtiTmp))
624 return false;
625
626 int const nXTmp = PtiTmp.nGetX();
627 int const nYTmp = PtiTmp.nGetY();
628
629 bool bPreExistingNotch;
630
631 // Get the existing notch apex elevation, if there is one
632 double dNotchApexElev = m_pRasterGrid->m_Cell[nXTmp][nYTmp].pGetLandform()->dGetCliffNotchApexElev();
633
634 if (bFPIsEqual(dNotchApexElev, DBL_NODATA, TOLERANCE))
635 {
636 // There is no existing notch apex elevation
637 bPreExistingNotch = false;
638
639 // So use the notch elevation for this timestep
641 }
642 else
643 bPreExistingNotch = true;
644
645 // So, can we create a notch here?
646 double const dSedTopElevNoTalus = m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetAllSedTopElevOmitTalus();
647 if (dNotchApexElev < dSedTopElevNoTalus)
648 {
649 // Yes we can potentially create a notch here
650#ifdef _DEBUG
651 if (bPreExistingNotch)
652 LogStream << m_ulIter << ": \tIncision of pre-existing inland cliff at [" << nXTmp << "][" << nYTmp << "] dNotchApexElev = " << dNotchApexElev << " dSedTopElevNoTalus = " << dSedTopElevNoTalus << endl;
653 else
654 LogStream << m_ulIter << ": \tCreation of new inland cliff at [" << nXTmp << "][" << nYTmp << "] dNotchApexElev = " << dNotchApexElev << " dSedTopElevNoTalus = " << dSedTopElevNoTalus << endl;
655#endif
656
657 // Set the cell to be an inland cliff
658 m_pRasterGrid->m_Cell[nXTmp][nYTmp].pGetLandform()->SetLFCategory(LF_CLIFF_INLAND);
659
660 // Set its apex elevation
661 m_pRasterGrid->m_Cell[nXTmp][nYTmp].pGetLandform()->SetCliffNotchApexElev(dNotchApexElev);
662
663 // We have some notch incision of this newly-created notch
664 double dWeight;
665 if (dWaveElev > dNotchApexElev)
666 // Should not happen: safety check
667 dWeight = 1;
668 else
669 // Assume a linear decrease in incision with distance downwards from notch apex
670 dWeight = 1 - ((dNotchApexElev - dWaveElev) / CLIFF_NOTCH_CUTOFF_DISTANCE);
671
672 // Calculate this-timestep cliff notch erosion (is a length in external CRS units). Note that only consolidated sediment can have a cliff notch
673 double const dNotchIncision = dWeight * dWaveEnergy / m_dCliffErosionResistance;
674
675 if (bPreExistingNotch)
676 m_pRasterGrid->m_Cell[nXTmp][nYTmp].pGetLandform()->AddToCliffNotchIncisionDepth(dNotchIncision);
677 else
678 m_pRasterGrid->m_Cell[nXTmp][nYTmp].pGetLandform()->SetCliffNotchIncisionDepth(dNotchIncision);
679
680 // And add to the cell's accumulated wave energy
681 m_pRasterGrid->m_Cell[nXTmp][nYTmp].pGetLandform()->AddToAccumWaveEnergy(dWaveEnergy * dWeight);
682
683#ifdef _DEBUG
684 double const dSedTopElevIncTalus = m_pRasterGrid->m_Cell[nXTmp][nYTmp].dGetAllSedTopElevIncTalus();
685
686 assert(dNotchApexElev < dSedTopElevNoTalus);
687
688 LogStream << m_ulIter << ": \tINLAND cliff created (or re-created) at [" << nXTmp << "][" << nYTmp << "] dNotchApexElev = " << dNotchApexElev << " dSedTopElevNoTalus = " << dSedTopElevNoTalus << " dSedTopElevIncTalus = " << dSedTopElevIncTalus << " dNotchIncision = " << dNotchIncision << " dNotchApexElev = " << dNotchApexElev << endl;
689#endif
690 bFound = true;
691 }
692
693 n++;
694 } while (! bFound);
695
696 // Should never get here
697 return bFound;
698}
699
700//===============================================================================================================================
702//===============================================================================================================================
703void CSimulation::DoCliffCollapseTalusDeposition(int const nCoast, CRWCliff const* pCliff, double const dSandFromCollapse, double const dCoarseFromCollapse, int const nNotchLayer)
704{
705 // Check: is there some sand- or coarse-sized sediment to deposit?
706 if ((dSandFromCollapse + dCoarseFromCollapse) < SED_ELEV_TOLERANCE)
707 return;
708
709 // Get the cliff cell's grid coords
710 int const nX = pCliff->pPtiGetCellMarkedAsCliff()->nGetX();
711 int const nY = pCliff->pPtiGetCellMarkedAsCliff()->nGetY();
712
713 // Get a pointer to the layer in which the notch was incised
714 CRWCellLayer* pLayer = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nNotchLayer);
715
716 // And get a pointer to the cell layer's talus object
717 CRWCellTalus* pTalus = pLayer->pGetOrCreateTalus();
718
719 if (dSandFromCollapse > 0)
720 {
721 // Add the sand-sized sediment from the collapse to the talus object for this layer
722 pTalus->AddSandDepth(dSandFromCollapse);
723 m_pRasterGrid->m_Cell[nX][nY].AddSandTalusDeposition(dSandFromCollapse);
724 }
725
726 if (dCoarseFromCollapse > 0)
727 {
728 // Add the coarse-sized sediment from the collapse to the talus object for this layer
729 pTalus->AddCoarseDepth(dCoarseFromCollapse);
730 m_pRasterGrid->m_Cell[nX][nY].AddCoarseTalusDeposition(dCoarseFromCollapse);
731 }
732
733 // And update the cell's sea depth
734 m_pRasterGrid->m_Cell[nX][nY].SetSeaDepth();
735
736#ifdef _DEBUG
737 LogStream << m_ulIter << ";\tcoast " << nCoast << " cliff collapse talus deposition on [" << nX << "][" << nY << "] dSandFromCollapse = " << dSandFromCollapse << " dCoarseFromCollapse = " << dCoarseFromCollapse << " sea depth = " << m_pRasterGrid->m_Cell[nX][nY].dGetSeaDepth() << endl;
738#endif
739}
740
741//===============================================================================================================================
743//===============================================================================================================================
745{
746 for (int nX = 0; nX < m_nXGridSize; nX++)
747 {
748 for (int nY = 0; nY < m_nYGridSize; nY++)
749 {
750 int const nLayers = m_pRasterGrid->m_Cell[nX][nY].nGetNumLayers();
751 for (int nLayer = 0; nLayer < nLayers; nLayer++)
752 {
753 // Is there talus on this cell layer?
754 CRWCellTalus* pTalus = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nLayer)->pGetTalus();
755
756 if (pTalus == NULL)
757 // No talus
758 continue;
759
760 // OK we have some talus which could be redistributed from this cell, if waves (inc runup) reach high enough
761 double const dThisTalusBottomElev = m_pRasterGrid->m_Cell[nX][nY].dGetAllSedTopElevOmitTalus();
762
763 // Find the nearest point on any coastline, so we can get the runup at this point
764 int nCoastClosest;
765 int const nCoastPointClosest = nFindClosestCoastPoint(nX, nY, nCoastClosest);
766 if (nCoastPointClosest == INT_NODATA)
768
769 // And get the wave runup at this point
770 double const dRunup = m_VCoast[nCoastClosest].dGetRunUp(nCoastPointClosest);
771
772 // Now calc the elevation to which waves reach
773 double const dWaveElev = m_dThisIterSWL + dRunup;
774
775 // Only move talus if waves (inc runup) reach above the bottom of the talus
776 if (dWaveElev < dThisTalusBottomElev)
777 {
778 // No talus moved
779#ifdef _DEBUG
780 LogStream << m_ulIter << ": \tNO talus moved from [" << nX << "][" << nY << "] since waves do not reach talus base: dWaveElev = " << dWaveElev << " dThisTalusBottomElev = " << dThisTalusBottomElev << endl;
781#endif
782 continue;
783 }
784
785 // OK we will move some talus
786 double const dThisTalusTopElev = m_pRasterGrid->m_Cell[nX][nY].dGetAllSedTopElevIncTalus();
787 double const dTalusDepth = dThisTalusTopElev - dThisTalusBottomElev;
788 double dWeight;
789
790 if (dWaveElev > dThisTalusTopElev)
791 // Should not happen: safety check
792 dWeight = 1;
793 else
794 // Assume a linear decrease in talus removal with distance downwards from talus top
795 dWeight = 1 - ((dThisTalusTopElev - dWaveElev) / dTalusDepth);
796
797 if (bFPIsEqual(dWeight, 0.0, TOLERANCE))
798 {
799#ifdef _DEBUG
800 LogStream << m_ulIter << ": \tNO talus moved from [" << nX << "][" << nY << "] dWeight = " << dWeight << endl;
801#endif
802 continue;
803 }
804
805#ifdef _DEBUG
806 LogStream << m_ulIter << ": \ttalus potentially moved from [" << nX << "][" << nY << "] dThisTalusBottomElev = " << dThisTalusBottomElev << " dThisTalusTopElev = " << dThisTalusTopElev << " dWeight = " << dWeight << endl;
807#endif
808
809 // Next, determine the cells to which talus will be moved. Find all surrounding cells with a top elevation (including talus) which is less than the top elevation (including talus) of this cell
810 double dAdjElev;
811 double dTotElevDiff = 0;
812 vector<double> VdAdjElevDiff;
813 vector<CGeom2DIPoint> VptAdj;
814
815 for (int nSearchDirection = NORTH; nSearchDirection <= NORTH_WEST; nSearchDirection++)
816 {
817 int nXAdj;
818 int nYAdj;
819
820 switch (nSearchDirection)
821 {
822 case NORTH:
823 nXAdj = nX - 1;
824 nYAdj = nY;
825
826 if (bIsWithinValidGrid(nXAdj, nYAdj))
827 {
828 dAdjElev = m_pRasterGrid->m_Cell[nXAdj][nYAdj].dGetAllSedTopElevIncTalus();
829 if (dAdjElev < dThisTalusTopElev)
830 {
831 VptAdj.push_back(CGeom2DIPoint(nXAdj, nYAdj));
832 VdAdjElevDiff.push_back(dAdjElev);
833 dTotElevDiff += dAdjElev;
834 }
835 }
836
837 break;
838
839 case NORTH_EAST:
840 nXAdj = nX;
841 nYAdj = nY - 1;
842
843 if (bIsWithinValidGrid(nXAdj, nYAdj))
844 {
845 dAdjElev = m_pRasterGrid->m_Cell[nXAdj][nYAdj].dGetAllSedTopElevIncTalus();
846 if (dAdjElev < dThisTalusTopElev)
847 {
848 VptAdj.push_back(CGeom2DIPoint(nXAdj, nYAdj));
849 VdAdjElevDiff.push_back(dAdjElev);
850 dTotElevDiff += dAdjElev;
851 }
852 }
853
854 break;
855
856 case EAST:
857 nXAdj = nX;
858 nYAdj = nY - 1;
859
860 if (bIsWithinValidGrid(nXAdj, nYAdj))
861 {
862 dAdjElev = m_pRasterGrid->m_Cell[nXAdj][nYAdj].dGetAllSedTopElevIncTalus();
863 if (dAdjElev < dThisTalusTopElev)
864 {
865 VptAdj.push_back(CGeom2DIPoint(nXAdj, nYAdj));
866 VdAdjElevDiff.push_back(dAdjElev);
867 dTotElevDiff += dAdjElev;
868 }
869 }
870
871 break;
872
873 case SOUTH_EAST:
874 nXAdj = nX + 1;
875 nYAdj = nY;
876
877 if (bIsWithinValidGrid(nXAdj, nYAdj))
878 {
879 dAdjElev = m_pRasterGrid->m_Cell[nXAdj][nYAdj].dGetAllSedTopElevIncTalus();
880 if (dAdjElev < dThisTalusTopElev)
881 {
882 VptAdj.push_back(CGeom2DIPoint(nXAdj, nYAdj));
883 VdAdjElevDiff.push_back(dAdjElev);
884 dTotElevDiff += dAdjElev;
885 }
886 }
887
888 break;
889
890 case SOUTH:
891 nXAdj = nX + 1;
892 nYAdj = nY;
893
894 if (bIsWithinValidGrid(nXAdj, nYAdj))
895 {
896 dAdjElev = m_pRasterGrid->m_Cell[nXAdj][nYAdj].dGetAllSedTopElevIncTalus();
897 if (dAdjElev < dThisTalusTopElev)
898 {
899 VptAdj.push_back(CGeom2DIPoint(nXAdj, nYAdj));
900 VdAdjElevDiff.push_back(dAdjElev);
901 dTotElevDiff += dAdjElev;
902 }
903 }
904
905 break;
906
907 case SOUTH_WEST:
908 nXAdj = nX + 1;
909 nYAdj = nY;
910
911 if (bIsWithinValidGrid(nXAdj, nYAdj))
912 {
913 dAdjElev = m_pRasterGrid->m_Cell[nXAdj][nYAdj].dGetAllSedTopElevIncTalus();
914 if (dAdjElev < dThisTalusTopElev)
915 {
916 VptAdj.push_back(CGeom2DIPoint(nXAdj, nYAdj));
917 VdAdjElevDiff.push_back(dAdjElev);
918 dTotElevDiff += dAdjElev;
919 }
920 }
921
922 break;
923
924 case WEST:
925 nXAdj = nX;
926 nYAdj = nY + 1;
927
928 if (bIsWithinValidGrid(nXAdj, nYAdj))
929 {
930 dAdjElev = m_pRasterGrid->m_Cell[nXAdj][nYAdj].dGetAllSedTopElevIncTalus();
931 if (dAdjElev < dThisTalusTopElev)
932 {
933 VptAdj.push_back(CGeom2DIPoint(nXAdj, nYAdj));
934 VdAdjElevDiff.push_back(dAdjElev);
935 dTotElevDiff += dAdjElev;
936 }
937 }
938
939 break;
940
941 case NORTH_WEST:
942 nXAdj = nX;
943 nYAdj = nY + 1;
944
945 if (bIsWithinValidGrid(nXAdj, nYAdj))
946 {
947 dAdjElev = m_pRasterGrid->m_Cell[nXAdj][nYAdj].dGetAllSedTopElevIncTalus();
948 if (dAdjElev < dThisTalusTopElev)
949 {
950 VptAdj.push_back(CGeom2DIPoint(nXAdj, nYAdj));
951 VdAdjElevDiff.push_back(dAdjElev);
952 dTotElevDiff += dAdjElev;
953 }
954 }
955
956 break;
957 }
958 }
959
960 int const nLower = static_cast<int>(VptAdj.size());
961 if (nLower == 0)
962 {
963 // None of the adjacent cells are lower
964#ifdef _DEBUG
965 LogStream << m_ulIter << ": \tNO talus moved from [" << nX << "][" << nY << "] since no adjacent cells are lower" << endl;
966#endif
967 continue;
968 }
969
970 // OK, at least one adjacent cell is lower. Move talus to each adjacent cell in proportion to the elevation difference
971 vector<double> VdPropToMove(nLower);
972 for (int n = 0; n < nLower; n++)
973 VdPropToMove[n] = VdAdjElevDiff[n] / dTotElevDiff;
974
975 // TODO Removal rate:
976 // * to be different for sand and coarse
977 // * to include talus erodibility
978 // * must depend on SWL and runup.
979 // Note we are ignoring subaerial processes.
980
981 double const dTalusSandOrig = pTalus->dGetSandDepth();
982 double const dTalusCoarseOrig = pTalus->dGetCoarseDepth();
983 double dTalusSandToMove = pTalus->dGetSandDepth();
984 double dTalusCoarseToMove = pTalus->dGetCoarseDepth();
985 double dTalusSandMoved = 0;
986 double dTalusCoarseMoved = 0;
987 double const dTalusErodibility = 0.3;
988 double const dSandRemovalRate = 1 * dTalusErodibility; // TEST external CRS units per hour e.g. metres depth per hour (since timestep is in hours)
989 double const dCoarseRemovalRate = 0.8 * dTalusErodibility; // TEST external CRS units per hour e.g. metres depth per hour (since timestep is in hours)
990
991 for (int n = 0; n < nLower; n++)
992 {
993 if (dTalusSandToMove > 0)
994 {
995 // We will deposit some talus sand onto the top layer of this adjacent cell
996 int const nXAdj = VptAdj[n].nGetX();
997 int const nYAdj = VptAdj[n].nGetY();
998
999 int const nTopLayer = m_pRasterGrid->m_Cell[nXAdj][nYAdj].nGetNumOfTopLayerAboveBasement();
1000 double const dSandOnAdj = m_pRasterGrid->m_Cell[nXAdj][nYAdj].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->dGetSandDepth();
1001
1002 double const dPotentialDepthToMove = pTalus->dGetSandDepth() * dWeight * dSandRemovalRate * VdPropToMove[n] * m_dTimeStep;
1003 double const dActualDepthToMove = tMin(dTalusSandToMove, dPotentialDepthToMove);
1004
1005 m_pRasterGrid->m_Cell[nXAdj][nYAdj].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->SetSandDepth(dSandOnAdj + dActualDepthToMove);
1006 dTalusSandToMove -= dActualDepthToMove;
1007 dTalusSandMoved += dActualDepthToMove;
1008
1009 assert(dTalusSandToMove >= 0.0);
1010
1011 // Set the changed-this-timestep switch re. the adjacent cell
1012 m_bUnconsChangedThisIter[nTopLayer] = true;
1013#ifdef _DEBUG
1014 LogStream << m_ulIter << ": \t" << std::scientific << dActualDepthToMove << std::fixed << " talus sand deposited at [" << nXAdj << "][" << nYAdj << "], talus sand still to deposit on [" << nX << "][" << nY << "] = " << std::scientific << dTalusSandToMove << " talus sand removed = " << dTalusSandMoved << std::fixed << endl;
1015#endif
1016 // TODO Update the adjacent cell's talus deposition, and total talus deposition, values
1017 // m_pRasterGrid->m_Cell[nX][nY].IncrBeachDeposition(dActualDepthToMove);
1018 }
1019
1020 if (dTalusCoarseToMove > 0)
1021 {
1022 // We will deposit some talus coarse onto the top layer of this adjacent cell
1023 int const nXAdj = VptAdj[n].nGetX();
1024 int const nYAdj = VptAdj[n].nGetY();
1025
1026 int const nTopLayer = m_pRasterGrid->m_Cell[nXAdj][nYAdj].nGetNumOfTopLayerAboveBasement();
1027 double const dCoarseOnAdj = m_pRasterGrid->m_Cell[nXAdj][nYAdj].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->dGetCoarseDepth();
1028
1029 double const dPotentialDepthToMove = pTalus->dGetCoarseDepth() * dWeight * dCoarseRemovalRate * VdPropToMove[n] * m_dTimeStep;
1030 double const dActualDepthToMove = tMin(dTalusCoarseToMove, dPotentialDepthToMove);
1031
1032 m_pRasterGrid->m_Cell[nXAdj][nYAdj].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->SetSandDepth(dCoarseOnAdj + dActualDepthToMove);
1033 dTalusCoarseToMove -= dActualDepthToMove;
1034 dTalusCoarseMoved += dActualDepthToMove;
1035
1036 assert(dTalusCoarseToMove >= 0.0);
1037#ifdef _DEBUG
1038 LogStream << m_ulIter << ": \t" << std::scientific << dActualDepthToMove << std::fixed << " talus coarse deposited at [" << nXAdj << "][" << nYAdj << "], talus coarse still to deposit on [" << nX << "][" << nY << "] = " << std::scientific << dTalusCoarseToMove << " talus coarse removed = " << dTalusCoarseMoved << std::fixed << endl;
1039#endif
1040
1041 // Set the changed-this-timestep switch re. the adjacent cell
1042 m_bUnconsChangedThisIter[nTopLayer] = true;
1043
1044 // TODO Update the adjacent cell's talus deposition, and total talus deposition, values
1045 // m_pRasterGrid->m_Cell[nX][nY].IncrBeachDeposition(dActualDepthToMove);
1046 }
1047 }
1048
1049 if (dTalusSandMoved > 0)
1050 {
1051 // For the source cell, update the sand talus value
1052 double const dTalusSandRemaining = tMax(dTalusSandOrig - dTalusSandMoved, 0.0);
1053
1054 pTalus->SetSandDepth(dTalusSandRemaining);
1055 }
1056
1057 if (dTalusCoarseMoved > 0)
1058 {
1059 // For this cell, update the cell layer's coarse talus value
1060 double const dTalusCoarseRemaining = tMax(dTalusCoarseOrig - dTalusCoarseMoved, 0.0);
1061
1062 pTalus->SetCoarseDepth(dTalusCoarseRemaining);
1063 }
1064
1065 // Has all the talus gone from this layer? If so, then delete it
1066 if (bFPIsEqual(pTalus->dGetSandDepth() + pTalus->dGetCoarseDepth(), 0.0, TOLERANCE))
1067 {
1068 LogStream << m_ulIter << ": \tpTalus->dGetSandDepth() + pTalus->dGetCoarseDepth() = " << std::scientific << pTalus->dGetSandDepth() + pTalus->dGetCoarseDepth() << std::fixed << " so deleting talus object at [" << nX << "][" << nY << "]" << endl;
1069 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nLayer)->DeleteTalus();
1070 }
1071
1072 // And update the cell's sea depth
1073 m_pRasterGrid->m_Cell[nX][nY].SetSeaDepth();
1074
1075#ifdef _DEBUG
1076 LogStream << m_ulIter << ": \ttalus moved from [" << nX << "][" << nY << "] sea depth = " << m_pRasterGrid->m_Cell[nX][nY].dGetSeaDepth() << endl;
1077#endif
1078 }
1079 }
1080 }
1081
1082 return RTN_OK;
1083}
Contains CGeom2DIPoint definitions.
Abstract class, used as a base class for landform objects on the coastline.
CGeom2DIPoint * pPtiGetCellMarkedAsCliff(void) const
Get the grid coordinates of the cell on which this cliff sits.
int nGetPointOnCoast(void) const
Get the point on the coast on which this coast landform sits.
int nGetLandFormCategory(void) const
Get the landform category.
void IncTotAccumWaveEnergy(double const)
Increment total accumulated wave energy.
Geometry class used to represent 2D point objects with integer coordinates.
Definition 2di_point.h:25
int nGetY(void) const
Returns the CGeom2DIPoint object's integer Y coordinate.
Definition 2di_point.cpp:51
int nGetX(void) const
Returns the CGeom2DIPoint object's integer X coordinate.
Definition 2di_point.cpp:45
Geometry class used for coast polygon objects.
void AddCliffCollapseToSuspensionFine(double const)
Add to the this-iteration total of unconsolidated fine sediment from cliff collapse which goes to sus...
void AddCliffCollapseErosionCoarse(double const)
Add to the this-iteration total of unconsolidated coarse sediment from cliff collapse on this polygon...
void AddCliffCollapseErosionFine(double const)
Add to the this-iteration total of unconsolidated fine sediment eroded from cliff collapse on this po...
void AddCliffCollapseErosionSand(double const)
Add to the this-iteration total of unconsolidated sand sediment from cliff collapse on this polygon,...
Real-world class used to represent the sediment layers associated with a cell object.
Definition cell_layer.h:28
CRWCellTalus * pGetOrCreateTalus(void)
Returns a pointer to the layer's talus object. If there is no talus object, then create one.
Real-world class used to represent the talus (unconsolidated sediment resulting from cliff collapse) ...
Definition cell_talus.h:25
double dGetCoarseDepth(void) const
Returns the coarse sediment depth equivalent for this talus object object.
double dGetSandDepth(void) const
Returns the sand sediment depth equivalent for this talus object.
void AddCoarseDepth(double const)
Adds coarse sediment (depth equivalent) to this talus object object's coarse sediment.
void SetSandDepth(double const)
Sets this talus object's sand sediment depth equivalent. Note no checks here to see if new equiv dept...
void SetCoarseDepth(double const)
Sets this talus object object's coarse sediment depth equivalent. Note no checks here to see if new e...
void AddSandDepth(double const)
Adds sand sediment (depth equivalent) to this talus object object's sand sediment.
Real-world class used to represent the 'cliff' category of coastal landform object.
Definition cliff.h:28
void SetNotchApexElev(double const)
Sets the elevation of the apex of the erosional notch (in external CRS units)
Definition cliff.cpp:70
bool bReadyToCollapse(double const) const
Returns true if the horizontal incision of the erosional notch exceeds the critical notch incision.
Definition cliff.cpp:99
void IncreaseNotchIncision(double const)
Increases the horizontal incision (in external CRS units) of the erosional notch, measured inland fro...
Definition cliff.cpp:88
void SetCliffCollapsed(void)
Flags the cliff as having collapsed.
Definition cliff.cpp:58
void SetNotchIncision(double const)
Sets the horizontal incision (in external CRS units) of the erosional notch, measured inland from the...
Definition cliff.cpp:76
double dGetNotchApexElev(void) const
Returns the elevation of the apex of the erosional notch (in external CRS units)
Definition cliff.cpp:64
int m_nLogFileDetail
The level of detail in the log file output. Can be LOG_FILE_LOW_DETAIL, LOG_FILE_MIDDLE_DETAIL,...
Definition simulation.h:588
int m_nNumTotCliffCollapse
The total number of cells with cliff collapse since the start of the simulation.
Definition simulation.h:537
double m_dThisIterSWL
The still water level for this timestep (this includes tidal changes and any long-term SWL change)
Definition simulation.h:726
CGeomRasterGrid * m_pRasterGrid
Pointer to the raster grid object.
int nMoveCliffTalusToUnconsolidated(void)
Move talus from previous cliff collapse to unconsolidated sediment.
int m_nXGridSize
The size of the grid in the x direction.
Definition simulation.h:477
ofstream LogStream
vector< CRWCoast > m_VCoast
The coastline objects.
int m_nNumThisIterCliffCollapse
The number of cells with cliff collapse this iteration.
Definition simulation.h:534
int m_nYGridSize
The size of the grid in the y direction.
Definition simulation.h:480
static CGeom2DIPoint PtiGetPerpendicular(CGeom2DIPoint const *, CGeom2DIPoint const *, double const, int const)
Returns a CGeom2DIPoint (grid CRS) which is the 'other' point of a two-point vector passing through P...
double m_dNotchIncisionAtCollapse
Notch overhang (i.e. length of horizontal incision at the apex elevation) to initiate collapse (m)
Definition simulation.h:927
void DoCliffCollapseTalusDeposition(int const, CRWCliff const *, double const, double const, int const)
Deposit the unconsolidated sediment from cliff collapse as talus on the cell on which collapse occurr...
double m_dThisIterCliffCollapseErosionCoarseUncons
This-iteration total of coarse unconsolidated sediment produced by cliff collapse (m^3)
Definition simulation.h:954
bool bCreateNotchInland(int const, int const, int const, int const, double const, double const)
If possible, creates an erosional notch further inland from a given coastline point.
double m_dThisIterCliffCollapseErosionFineUncons
This-iteration total of fine unconsolidated sediment produced by cliff collapse (m^3)
Definition simulation.h:948
double m_dThisIterNewNotchApexElev
Elevation (m) of the apex of any cliff notches created during this iteration.
Definition simulation.h:930
double m_dCliffErosionResistance
Resistance of cliff to notch erosion.
Definition simulation.h:924
double dGridCentroidYToExtCRSY(int const) const
Given the integer Y-axis ordinate of a cell in the raster grid CRS, returns the external CRS Y-axis o...
Definition gis_utils.cpp:74
int nDoCliffCollapse(int const, CRWCliff *, double &, double &, double &, int &, double &, double &)
Simulates cliff collapse on a single cell. Collapse happens when when a notch which is incised into t...
bool bIncreaseCliffNotchIncision(int const, int const, int const, CRWCliff *, double const)
Increase the incision (if any) of a cliff notch, assuming a linear decrease in incision with distance...
double m_dThisIterUnconsCoarseCliffDeposition
This-iteration total of coarse unconsolidated sediment deposited due to cliff collapse (m^3)
Definition simulation.h:969
double m_dThisIterUnconsSandCliffDeposition
This-iteration total of sand unconsolidated sediment deposited due to cliff collapse (m^3)
Definition simulation.h:966
double m_dThisIterFineSedimentToSuspension
Total fine unconsolidated sediment in suspension for this iteration (depth in m)
Definition simulation.h:879
int nDoAllWaveEnergyToCoastLandforms(void)
double m_dNotchApexAboveMHW
Distance of notch base below SWL (m)
Definition simulation.h:933
vector< bool > m_bUnconsChangedThisIter
One element per layer: has the consolidated sediment of this layer been changed during this iteration...
double m_dThisIterCliffCollapseErosionFineCons
This-iteration total of fine consolidated sediment produced by cliff collapse (m^3)
Definition simulation.h:957
double m_dTimeStep
The length of an iteration (a timestep) in hours.
Definition simulation.h:690
double m_dCellArea
Area of a cell (in external CRS units)
Definition simulation.h:675
bool bIsWithinValidGrid(int const, int const) const
Checks whether the supplied point (an x-y pair, in the grid CRS) is within the raster grid,...
int nFindClosestCoastPoint(int const, int const, int &)
Finds the number of the closest point on any coastline to a given point, or INT_NODATA in case of err...
double m_dThisIterCliffCollapseErosionSandCons
This-iteration total of sand consolidated sediment produced by cliff collapse (m^3)
Definition simulation.h:960
double m_dThisIterCliffCollapseErosionSandUncons
This-iteration total of sand unconsolidated sediment produced by cliff collapse (m^3)
Definition simulation.h:951
unsigned long m_ulIter
The number of the current iteration (time step)
Definition simulation.h:609
double dGridCentroidXToExtCRSX(int const) const
Given the integer X-axis ordinate of a cell in the raster grid CRS, returns the external CRS X-axis o...
Definition gis_utils.cpp:65
double m_dCellSide
Length of a cell side (in external CRS units)
Definition simulation.h:672
double m_dThisIterCliffCollapseErosionCoarseCons
This-iteration total of coarse consolidated sediment produced by cliff collapse (m^3)
Definition simulation.h:963
vector< bool > m_bConsChangedThisIter
One element per layer: has the consolidated sediment of this layer been changed during this iteration...
Contains CRWCliff definitions.
This file contains global definitions for CoastalME.
int const INT_NODATA
Definition cme.h:380
T tMin(T a, T b)
Definition cme.h:1175
double const TOLERANCE
Definition cme.h:725
int const RTN_ERR_NO_TOP_LAYER
Definition cme.h:629
int const ELEV_ABOVE_SEDIMENT_TOP
Definition cme.h:671
int const SOUTH_EAST
Definition cme.h:402
int const LOG_FILE_MIDDLE_DETAIL
Definition cme.h:393
bool bFPIsEqual(const T d1, const T d2, const T dEpsilon)
Definition cme.h:1213
int const LF_CLIFF_INLAND
Definition cme.h:438
int const SOUTH_WEST
Definition cme.h:404
T tMax(T a, T b)
Definition cme.h:1162
int const NORTH_WEST
Definition cme.h:406
int const RTN_ERR_CLIFF_NOT_IN_POLYGON
Definition cme.h:652
string const WARN
Definition cme.h:806
int const SOUTH
Definition cme.h:403
double const CLIFF_NOTCH_CUTOFF_DISTANCE
Definition cme.h:735
int const LOG_FILE_ALL
Definition cme.h:395
int const EAST
Definition cme.h:401
int const RTN_OK
Definition cme.h:585
int const LF_CLIFF_ON_COASTLINE
Definition cme.h:437
int const NORTH_EAST
Definition cme.h:400
double const DBL_NODATA
Definition cme.h:736
int const RTN_ERR_CLIFF_NOTCH
Definition cme.h:621
int const ELEV_IN_BASEMENT
Definition cme.h:670
double const SED_ELEV_TOLERANCE
Definition cme.h:726
int const RTN_ERR_CLIFF_TALUS_TO_UNCONS
Definition cme.h:666
int const NORTH
Definition cme.h:399
int const WEST
Definition cme.h:405
Contains CACoastLandform definitions.
Contains CSimulation definitions.