MindTheGap-ERC / CarboKitten.jl

Julia implementation of carbonate platform model
https://mindthegap-erc.github.io/CarboKitten.jl/
GNU General Public License v3.0
1 stars 0 forks source link

Modeling subaerial erosion in kyr scale #16

Open xyl96 opened 10 months ago

xyl96 commented 10 months ago

Aims

The sub-aerial exposure of the carbonate platform could remove significant amount of strata previously deposited. For example, in Natih formation, ancient channels are observed (Grélaud et al., 2016). More recently, 35Cl data is deployed on recently exposed carbonate terrace and demonstrated that the denudations are indeed happening and under the situation of monsoon tropical climate, the rate is estimated to be around 20mm/kyr (Chauveau et al., 2021) However, most carbonate architecture models did not take this into account due to the complexity of denudation.

Herein, we have two ways to quantify the sub-aerial denudation: 1) empirical relationship derived from 35Cl data and 2) conducting metamodels from landscape evolution models (e.g., WEPP, GOLIME, CAESAR, etc.). The first approach solely considers the remove of materials in a certain time period while the second approach could consder both removal of material and redeposition. We could compare the outcomes of the two different approaches with the original CarbonKitten (i.e., no denudation) and/or with the constant denudation rates. Ideally, we could utilise mathmatical tools to quantitatively tell the differences among them.

Empirical relationship derived from 35Cl data

Research based on the karst region and carbonate platform terrace suggested that the denudation rates are mainly controlled by precipitation and slopes, although the debates about which factor is more important is still ongoing (Ryb et al., 2014, Thomas et al., 2018). In general, the precipitation mainly controls the chemical dissolution while the slopes mainly controls the physical ersions. In addition, the type of carbonates may also play an important role (Kirklec et al., 2022), but given this feature is not studied well we will ditch it for now. Yang et al., 2020 have compiled the denudation rates (mm/kyr) along with precipitation and slopes and therefore we can use their data to create a function relates denudation rates (mm/kyr) to precipitation and slopes (their data is here). This is an empirical relationship and have a relatively large uncertainty in terms of fitting.

image

Fig 1. The relationship between MAP (mean precipitation per year, mm/y) and denudation rates (mm/ky)

image

Fig 2. The relationship between the curvature (i.e., slope or height) and the denudation rates (mm/ky)

We can see that both the slope and precipitation could increase the denudation rates, and reaches a 'steady state' after a certain point.

In terms of impletementation, we could simply feed the relationship here with kyr-averaged precipitation data, and the slope data. However, this emperical method may have issues in estimating sediment traps.

Forward modelling

This approach consists of two parts:

and the total denudation = chemical weathering + physical erosion.

The chemical weathering part has been documented in issue #5. While the physical erosion is listed below:

The presence of fluvial systems are debatable in carbonate platfrom tops due to highly permeable features of carbonate, but local redistribution of detached sediments by water are typically found (e.g., in Bahamas). That is to say, the platforms have limited long-range transportations for eroded carbonate sediments, but may still have transportation of materials to the neighboring place. Herein, I propose we could couple the cellular models proposed by Tucker et al., 1998 and van de Wiel et al., 2006.

In a nutshell, we first calculate how much material would be eroded in a cell based on the following equation:

image

Q is surficial run-off, S is gradient, m=1/3, n=2/3. The hb is the height. Q = (P-I) A, where P is precipitation, I is infiltration, and A is the area of a catchment. We can simply just apply an infiltration factor to P, and Q = P (1 - infiltration_factor), where this factor varies between 0 and 1. As we do not consider the fluvial systems, A = the area of one cell in this context. That is to say, only the run-off produced within the cell will be considered and by doing so we can avoid creating fluvial systems as we only consider the short-range erosion. As chemical weathering part suggested, The infiltration is difficult to estimate, worthy to have a discussion. Q is the max discharge (A P). Therefore Q/Q* is a rate depicting 'the fullness' of the run-off. kv = 0.23m/y in the paper, and the 'resistant rocks' are 0.023m/y. This is a very debatable constant! Worthy have a discussion. S is slope and the original paper defines it as the tan of slope angle Howard 1994. Herein, as only short-ranged water-based transportation is considered, we could only consider the slopes between two neigboring cells.

This equation comes from incision of fluivial systems into bedrocks. Thus, compared to equations desgined to describe eroding 'loose sediments', this is more appropriate to be used in carbonate platform where carbonate sediments are lithified quickly. Applying this equation does not mean we try to model the evolution of fluvial systems, instead, we only use this equation to calculate how much material is removed in a single cell.

Can we apply a equation designed to model river erosion in to a single cell? I think we can do here. This is because this equation derives from equations descrbing detaching rocks bu shear stress, and this equation has nothing to do with river 'morphology':

image

where the shear stress $\tau$ depends on Q (discharge) and slope (S). This makes sense: higher discharge and higher slope would allow higher water velocity and higher erosive ability. However, it is still questionable whether the run-off generated from a single cell could produce enough shear stress to detach lithified carbonate rocks.

Such short-ranged physical erosion is very similar to soil creeping and "diffusion' Kaufman et al., 2001:

image

where kd is a constant descrbing diffusivity. The diffusion equation is developed based on regolith mass-wasting and it is unclear whether we can apply it to bedrock-like carbonate rocks. However, we could still have a try and compare them.

The deposition of the eroded material is distributed to the neighboring 8 cells according to the differences in height van de Wiel et al., 2006.

image

Herein, the Vi is the amount of eroded material waiting to be transported to the neighboring lower cells = $\Delta$ hb * A (cell area). Vi,j is the amount of material that are recieved in each cell.

So We would run this simple model with time step of second (or if easier, in day?, worthy have a discussion), and utilise the meterologic data of Bahamas, and run it for 1kyr. The workflow would be in the follwoing pic: image Note that whether the coupling of chemical dissolution should be in the loop or not should be discussed. The initial topography could be randomly generated and/or dervied from CarboKitten.

Ideally, we could get a relationship between the topography and denudation rates in kyr scale: image Then we can apply this 'relationship' into the CarboKitten. Just like the production function, it will be a function of $\Delta$ h.

Alternatively, we can get a rule similar to CellularAutomata: we consider 9 cells (or 16 cells), and we derive a rule of the amount of erosion and deposition from the model described above: for a single cell, calculating how many cells are topographically higher than the target cell and how many are lower. We can assigned a weighted contribution factor based on the $\Delta$ h of the surrounding cells (differences in height), and caculate whether this cell would receive or erode (see the following diagram). Such weighted contribution factor is the cellular automata rule we are trying to derive from the model described above. This is slightly different from the cellular automata we used for the production, as the cellular automata rules we considered herein not only think about how many cells are neighboring, but also to what extent they are different. image For example, each cell has a vector to describe its topograhy data. the first number is it's initial topography at a certain time, the second is the eroded amount in one time step and the third is the deposted amount in one time step. The amount of materials being eroded and deposited depends on the topographic maps. By doing so, we not only extend the timescale from landscape evolution model from seconds or days to 1kyr, but also significantly reduce the calculations within the CarboKitten model by replacing the process-based model with simpler few math equations or Cellular automata rules.

Note this means that this relationship only would work under a certain climate conditions, for example, in this case we use Tropical Monson scenario (Bahamas).

What we gonna do now

We can start with the emperical data and try to use the emperical relationship to quantify the erosion rates based on slope and precipitation.

We can also start quantifying the chemical dissolution of rocks (#5).

We can begin to write codes of calculating gradient matrix.

EmiliaJarochowska commented 10 months ago

Thanks for this, @xyl96. Is it correct that it is the updated version after our discussion on Tuesday?

In the first equations: What is $k_v$? Over what area is $S$ calculated?

Can $P - I$ be simplified? By keeping both we introduce two new unknowns into the model. $P$ may seem easy to obtain, but will be hard to get for 1 ky timescales and I is even less constrained: can be near 0 to 100% $P$ depending on rock type. Instead of varying these two parameters we could just use a parameter that corresponds to the surface discharge?

This could have the advantage that if $P$ and I don’t vary between cells so they don’t need to be calculated per each cell?

However, it is still questionable whether the run-off generated from a single cell could produce enough shear stress to detach lithified carbonate rocks.

What are the critical values of shear stress to disintegrate lithified carbonates? Perhaps we will never reach them so we won’t need to worry about the second part?

Note that whether the coupling of chemical dissolution should be in the loop or not should be discussed.

Whatever is disintegrated will be subject to chemical erosion so I think it should be coupled. But another question is how this relates to chemical erosion of lithified rock: should it be applied first, and then mechanical erosion?

I cannot tell if the CA approach is not too complicated. The one you described first, based on discharge, seems simple and probably sufficient. The CA would add quite some extra computations. What did @jhidding say about that?

What we gonna do now

There are some things in this specification that you left open. So perhaps we should tie these up before you start coding?

  1. "The infiltration is difficult to estimate, worthy to have a discussion" - See my comment above, is it possible to simplify $P - I$?
  2. "kv = 0.23m/y in the paper, and the 'resistant rocks' are 0.023m/y. This is a very debatable constant! Worthy have a discussion." This is a parameter, not part of the algorithm, so I believe we can always modify it for a given run, so the discussion is not needed for the implementation.
  3. "We would run this simple model with time step of second (or if easier, in day?, worthy have a discussion)" I didn't get whether you want to apply it only if lithified rock is disintegrated into transportable sediment, because then maybe we never reach that stage?

We can also start quantifying the chemical dissolution of rocks (https://github.com/MindTheGap-ERC/CarboKitten.jl/issues/5).

Unless @jhidding disagrees, this is perhaps what you can indeed start while we clarify the points above?

NiklasHohmann commented 10 months ago
  1. @xyl96 you use the equation

    image

    to calculate erosion in a single cell. All cells are eroded simultaneously, so this will lead to a large coupled system of ODEs. Do you plan to (1) solve the whole system or (2) ignore the spatial dependency and only determine gradients every time step (see also comment below)? The consideration of spatial dependency also holds for redistribution of sediments via

    image

    How is this resolved?

  2. Because we use a rectangular grid, transport distances between neighboring cells differ (e.g., diagonal neighbors). Is a correction for this necessary?

  3. Regarding time steps: You mention time steps of seconds or days, and kyrs. Can you specify which parts of the model you would solve in continuous time (e.g. solve as ODE on whatever time steps the solver chooses), and which ones will be run on the time increments selected in CarboKitten?

xyl96 commented 10 months ago

Thanks for this, @xyl96. Is it correct that it is the updated version after our discussion on Tuesday?

In the first equations: What is kv? Over what area is S calculated?

Can P−I be simplified? By keeping both we introduce two new unknowns into the model. P may seem easy to obtain, but will be hard to get for 1 ky timescales and I is even less constrained: can be near 0 to 100% P depending on rock type. Instead of varying these two parameters we could just use a parameter that corresponds to the surface discharge? > Yes that what I thought as well - applying an infiltration factor on P. Updated! This could have the advantage that if P and I don’t vary between cells so they don’t need to be calculated per each cell?

However, it is still questionable whether the run-off generated from a single cell could produce enough shear stress to detach lithified carbonate rocks.

What are the critical values of shear stress to disintegrate lithified carbonates? Perhaps we will never reach them so we won’t need to worry about the second part? > OK!

Note that whether the coupling of chemical dissolution should be in the loop or not should be discussed.

Whatever is disintegrated will be subject to chemical erosion so I think it should be coupled. But another question is how this relates to chemical erosion of lithified rock: should it be applied first, and then mechanical erosion? > If we do not consider the relithification, I think first mechanical erosion, because chemical ones depends on the physical ones I cannot tell if the CA approach is not too complicated. The one you described first, based on discharge, seems simple and probably sufficient. The CA would add quite some extra computations. What did @jhidding say about that? > OK, I remebered on last Tuesday Johan also suggested that.

What we gonna do now

There are some things in this specification that you left open. So perhaps we should tie these up before you start coding?

  1. "The infiltration is difficult to estimate, worthy to have a discussion" - See my comment above, is it possible to simplify P−I?
  2. "kv = 0.23m/y in the paper, and the 'resistant rocks' are 0.023m/y. This is a very debatable constant! Worthy have a discussion." This is a parameter, not part of the algorithm, so I believe we can always modify it for a given run, so the discussion is not needed for the implementation.
  3. "We would run this simple model with time step of second (or if easier, in day?, worthy have a discussion)" I didn't get whether you want to apply it only if lithified rock is disintegrated into transportable sediment, because then maybe we never reach that stage? I am sorry, I think I did not fully get your question here. What I original mean is the original GOLIME have a short time steps, and if we use their equations we may also need to apply short-timescales

    We can also start quantifying the chemical dissolution of rocks (#5).

Unless @jhidding disagrees, this is perhaps what you can indeed start while we clarify the points above? Alright!

EmiliaJarochowska commented 10 months ago

Last week we discussed some of the questions in person. Currently the specification contains several different variants ("need to discuss", "we could do this" etc) so it is still not clear what decided to do, @xyl96? Is there any input from us that you are missing?

Above, you agreed that the approach based on CA is probably not necessary. But then you still discuss the time scale to run the model by Van de Wiel et al. (2007), or so it follows from the context:

So We would run this simple model with time step of second

This follows just after you describe CAESAR, but in your answer later you seem to mean GOLIME in this context:

I am sorry, I think I did not fully get your question here. What I original mean is the original GOLIME have a short time steps, and if we use their equations we may also need to apply short-timescales

Partly this might be the result of complicated formatting of your answer. If you use the citations correctly, you should be able to cite several levels while keeping it clear what the new text is.

May I suggest the following:

  1. Please do not edit the original description, because this will make following the discussion harder. Instead update the specification based on what we discussed, cutting out the things we decided to discard. You can mention them, but then clearly as not implemented here. Maybe we come back to them later, after this initial task is defined and completed.
  2. The text you wrote is a good literature review. It could probably be pasted to an article introduction. But for the purposes of what you are going to do, we do not need a list of possible variants, but just the one that you think is the best. Or if you need input on making a decision, please indicate this.
  3. @NiklasHohmann raised some important points, I am looking forward to hearing what you think about his questions.
xyl96 commented 10 months ago

Updated: writing it as a new comment is for cleaness and better visualization

Aims

The sub-aerial exposure of the carbonate platform could remove significant amount of strata previously deposited. For example, in Natih formation, ancient channels are observed (Grélaud et al., 2016). More recently, 35Cl data is deployed on recently exposed carbonate terrace and demonstrated that the denudations are indeed happening and under the situation of monsoon tropical climate, the rate is estimated to be around 20mm/kyr (Chauveau et al., 2021). That is to say, in each glacial cycle, ~2m thickness of carbonates would be removed (if assuming the cycle is 100kyr). However, most carbonate architecture models did not take this into account due to the complexity of denudation.

Herein, we have two ways to quantify the sub-aerial denudation: 1) empirical relationship derived from 35Cl data and 2) conducting metamodels from landscape evolution models (e.g., WEPP, GOLIME, CAESAR, etc.). The first approach solely considers the remove of materials in a certain time period while the second approach could consder both removal of material and redeposition. In addition, the observed data is also the dataset to constrain and calibrate the second approach. We could compare the outcomes of the two different approaches with the original CarbonKitten (i.e., no denudation) and/or with the constant denudation rates. Ideally, we could utilise mathmatical tools to quantitatively tell the differences among them.

Empirical relationship derived from 35Cl data

Research based on the karst region and carbonate platform terrace suggested that the denudation rates are mainly controlled by precipitation and slopes, although the debates about which factor is more important is still ongoing (Ryb et al., 2014, Thomas et al., 2018). In general, the precipitation mainly controls the chemical dissolution while the slopes mainly controls the physical ersions. In addition, the type of carbonates may also play an important role (Kirklec et al., 2022), but given this feature is not studied well we will ditch it for now. Yang et al., 2020 have compiled the denudation rates (mm/kyr) along with precipitation and slopes and therefore we can use their data to create a function relates denudation rates (mm/kyr) to precipitation and slopes (their data is here). This is an empirical relationship and have a relatively large uncertainty in terms of fitting.

image

Fig 1. The relationship between MAP (mean precipitation per year, mm/y) and denudation rates (mm/ky)

image

Fig 2. The relationship between the curvature (i.e., slope or height) and the denudation rates (mm/ky)

We can see that both the slope and precipitation could increase the denudation rates, and reaches a 'steady state' after a certain point.

In terms of impletementation, I used the function form of D = P * S, where D means denudation rates, P means effects of precipitation while S means effects of Slope. By doing so, we can consider both effects. Such formula structure is similar to RUSLE model, a widely used LEM (e.g., (Thapa et al., 2020)[https://environmentalsystemsresearch.springeropen.com/articles/10.1186/s40068-020-00177-2]). In function 'P', I used linear piece-wise function, as the Yang et al., 2020 paper suggested. For function 'S', I also used a piece-wise function as the figure shows: when slope is low, it's linear and when slope is high it's log.

Forward modelling

This approach consists of two parts:

and the total denudation = chemical weathering + physical erosion.

The chemical weathering part has been documented in issue #5. While the physical erosion is listed below:

Herein, I propose we could couple the cellular models proposed by Tucker et al., 1998 and van de Wiel et al., 2006.

In a nutshell, we first calculate how much material would be eroded in a cell based on the following equation (Tucker et al., 1998):

image

Q is surficial run-off, S is gradient, m=1/3, n=2/3. The hb is the height. Q = (P-I) A, where P is precipitation, I is infiltration, and A is the area of a catchment. We can simply just apply an infiltration factor to P, and Q = P (1 - infiltration_factor), where this factor varies between 0 and 1. As we do not consider the fluvial systems, A = the area of one cell in this context. That is to say, only the run-off produced within the cell will be considered and by doing so we can avoid creating fluvial systems as we only consider the short-range erosion. Q is the max discharge (A P). Therefore Q/Q* is a rate depicting 'the fullness' of the run-off. kv = 0.23m/y (at least for now). S is slope and the original paper defines it as the tan of slope angle Howard 1994 (more details please see below). Herein, as only short-ranged water-based transportation is considered.

This equation comes from incision of fluivial systems into bedrocks. Thus, compared to equations desgined to describe eroding 'loose sediments', this is more appropriate to be used in carbonate platform where carbonate sediments are lithified quickly. Applying this equation does not mean we try to model the evolution of fluvial systems, instead, we only use this equation to calculate how much material is removed in a single cell.

Can we apply a equation designed to model river erosion in to a single cell? I think we can do here. This is because this equation derives from equations descrbing detaching rocks bu shear stress, and this equation has nothing to do with river 'morphology':

image

where the shear stress $\tau$ depends on Q (discharge) and slope (S). This makes sense: higher discharge and higher slope would allow higher water velocity and higher erosive ability.

Such short-ranged physical erosion is very similar to soil creeping and "diffusion' Kaufman et al., 2001:

image

where kd is a constant descrbing diffusivity. The diffusion equation is developed based on regolith mass-wasting and it is unclear whether we can apply it to bedrock-like carbonate rocks. I will also try this equation and do some comparisons after finishing the one shown above.

The deposition of the eroded material is distributed to the neighboring 8 cells according to the differences in height van de Wiel et al., 2006.

image

Herein, the Vi is the amount of eroded material waiting to be transported to the neighboring lower cells = $\Delta$ hb * A (cell area). Vi,j is the amount of material that are recieved in each cell.

So We would run this simple model with time step the same as CarboKitten does. I would try 1) the default parameters in the literature first, 2) fitting the observed data in observed Cl isotopes and get parameters like kv (erodbility).

Implementation

The workflow would be in the follwoing pic: image The chemical dissolution is designed to be impleted after the physical erosion in each timestep. Each timestep = the timestep of CarboKitten.

Starting with one cell and its 8 neighbors, we first find out how many cells are lower than this cell. For example, if the targeting cell is the highest among the 9 cells, it suggests that this cell would be eroded and the eroded materials would be redistributed to the 8 neighboring cells. The composite slope would be calculated by weighting the slopes from 8 directions. Herein, I will give the diagonal cells weighting factor of 1, and direct neighboring cells of 2 ((reference)[https://pro.arcgis.com/en/pro-app/latest/tool-reference/spatial-analyst/how-slope-works.htm]. Therefore, in total we have 12 weights (42 + 41). We will use the weighted-averaged slope (tan of degrees) to calculate the eroded material. For now, we do not consider the global spatial heterogeneities, we only consider the cell and the neighboring 8 cells: this saves computational sources.

What we gonna do now

We can start with the emperical data and try to use the emperical relationship to quantify the erosion rates based on slope and precipitation. (doing right now)

We can also start quantifying the chemical dissolution of rocks (#5). (doing right now)

We can begin to write codes of calculating gradient matrix.

NiklasHohmann commented 10 months ago

Thanks for the update! Some comments

xyl96 commented 10 months ago

Thanks so much for these comments. I have updated it again.

  • RE volume of removed material: 20 mm/ kyr * 100 kyr = 2000 mm = 2 m removed.

Thanks for this!! Yeah my bad.

  • RE empirical relationship:

    1. You make a good point that slope and precipitation are driving factors of denudation. But what are the functional relationships involved? What functions will you use to model this relationship? This is the key part of the specification that must be as clear as possible. "simply feed the relationship here with kyr-averaged precipitation data, and the slope data" is not specific enough, there should be an equation on how input parameters (precipitation, slope) determine outputs (denudation rate)

Updated. Briefly, using piece-wise function. It's in the repositary Non-CarboKitten (which means everything else except the real hardcore CarboKitten stuff).

  1. The figures show denudation as a function of precipitation and of gradient independently. How you plan to identify how denudation rates change when both parameters are varied together?

Yeah, that's why I have to do regression on two parameters at the same time. The original dataset itsef have both slope and precipitation data, and this enables us to do so.

  • RE soil diffusion & creeping: Is there any indication that this is a relevant process on carbonate platforms?

I don't think so. But what I mean is we can vary the constants (diffussivity) in the diffusion eqn and see how it looks like (More like a trial). Diffusion just a way to smooth out the surface and from this point of view I think worthy to have a try.

  • RE "fitting the observed data in Cl isotopes and get parameters.":

    1. In my understanding you proposed two potential models: One based on the empirical relationship from Cl isotopes, one based on landscape evolution model. Do I understand correctly here that you want to calibrate the landscape evolution model with Cl isotope data?

Yes you are 100% right. But for the forward modeling we still need data to calibrate. The only data we had for 1kyr scale is the Cl isotope data (so far).

  1. Which parameters are you referring to here?

Updated. Bascially the constants in the equations.

xyl96 commented 7 months ago

Add tests for denudation ASAP.

EmiliaJarochowska commented 6 months ago

Thanks, @xyl96, for keeping tabs on it. Could you please make this a new issue? Please don't assign anyone yet, I would like to discuss it with @jhidding and @HannoSpreeuw in our meeting in the coming week and then we can assign the right person(s).