METHOD FOR MODELLING TRANSPORT AND DEPOSITION OF SEDIMENTS AND THE EVOLUTION OF A SEDIMENTARY BASIN

Abstract

A computer-implemented methods of modelling transport of a particle induced by water currents in an immersed area, and of modelling sedimentary deposition within the immersed area where particles represent sediments is described. The immersed area comprises a plurality of cells associated with respective water depths, and the particle is located in a cell and represents a quantity of sediments of determined granulometry and sediment type. The method includes determining a direction and velocity of at least one water current occurring within the immersed area, determining, from the water current, a direction and intensity of a shear stress induced by the water current on a particle, and determining a transport of the particle, from the determined direction and intensity of the shear stress induced on the particle, the granulometry and sediment type of the particle.

Description

BACKGROUND

Technical Field

The disclosure relates to a computer-implemented method for modelling water current-induced transport of particles of sediments in an immersed area, and to a computer implemented method for modelling the evolution of a sedimentary basin including transport and deposition of particles of sediments.

Description of the Related Art

Forward stratigraphic modelling is already known for modelling the evolution of sedimentary basins. In this type of modelling, an area is defined as a geological gridded model, and the modelling comprises superposing layers on the gridded model, each layer corresponding to a predetermined period of time and having a thickness which depends on an amount of material brought or created at a defined location during the period of time.

An example of forward stratigraphic modelling is for instance the DionisosFlow™ numerical stratigraphic model developed by IFP Energies Nouvelles, which allows reconstructing the stratigraphic architecture of sedimentary basins at a regional scale, by modelling basin deformation, clastic and carbonate supplies and sediment transport in continental and marine environment.

This model is used to simulate areas of regional scale, i.e., having dimensions of about tens to hundreds of kms side length², and on very long time scales. Each time layer simulated in this model is of at least 1000 years, up to 10,000 years, in order to model phenomena occurring on durations of at 100,000 years to several million years.

In order to be able to compute phenomena on such large geographical and time scales, the algorithms implemented by this model for simulating the transport of particles are only based on advection and diffusion. This implies that the geological phenomena which are simulated in forward simulators like DionisosFlow™ are necessarily continuous and homogeneous, and hence they cannot render local heterogeneities, which limits their potential for instance to precisely model the formation of oil reservoirs.

On the other hand, some models also exist to simulate phenomena which are more localized both in time and space, such as snow avalanches, sedimentation dynamics within rivers and deltas, etc. These models are based on equations modelling the motion of fluids such as Navier-Stokes equations. These models require a high computational load as compared to the time scale of the modelled phenomena, since the modelling of a phenomena requires a computational time of 1.5 times the duration of the modelled phenomena. Hence these times of models are not suitable for modelling the formation of reservoirs.

It is also known from WO 2020/229863, WO 2020/229862, WO 2020/229866 and WO 2020/229865 methods for modelling the evolution of sedimentary basins by simulating current-induced particle transport. These methods are forward stratigraphic modelling methods enabling to model the deposition of time-layer of sediments corresponding to an amount of time that is less than the duration of a time-layer of DionisosFlow™, for instance between a few 100 years and a few 10,000 years. The methods disclosed in these documents allow modelling the transport of particles based on current transport and takes into account the possible interaction between several phenomena or several types of current transport.

BRIEF SUMMARY

An aim of the disclosure is to provide a method for modelling the transport of particles and the evolution of sedimentary basins that is more heterogeneous and realistic than the prior art. Accordingly, a sedimentary basin output by the method can be used to extract a log of sediments that can be more accurately compared to an observed well.

Accordingly, a computer-implemented method of modelling transport of a particle induced by water currents in an immersed area is disclosed, wherein the immersed area comprises a plurality of cells associated with respective water depths, and the particle is located in a cell and represents a quantity of sediments of determined granulometry and sediment type, the method comprising:

• • determining a direction and velocity of at least one water current occurring within the immersed area,
  • determining, from the water current, a direction and intensity of a shear stress induced by the water current on a particle, and,
  • determining a transport of the particle, from the determined direction and intensity of the shear stress induced on the particle, the granulometry and sediment type of the particle.

In embodiments, the method further comprises a preliminary step of defining three water layers corresponding to respective water depths ranges extending between the surface and the bottom of the immersed area, comprising:

• • a bottom layer, located at water bottom,
  • a plume layer, located at water surface and,
  • a subsurface layer, extending between the bottom layer and the plume layer, and the method comprises determining the direction and intensity of a shear stress in at least one water layer.

In embodiments, the direction and intensity of a shear stress in a water layer is determined based on a velocity profile according to depth of the at least one water current in the considered water layer.

In embodiments, the method further comprises determining a direction and intensity of a shear stress induced by the gravity on particles located in the bottom layer.

In embodiments, the method the shear stress induced by the gravity on particles is determined such that:

• • the shear stress induced by gravity on the particle is null if the topographic slope of the water bottom is null,
  • the value of the shear stress induced by gravity on the particle is equal to a deposition shear stress threshold when the current velocity in the subsurface layer is null and the topographic slope of the water bottom is superior or equal to an avalanche angle determined for the particle, and
  • the direction of the shear stress induced by gravity on the particle is parallel to the direction of the downward topographic slope of the water bottom.

In embodiments, the particle is located within one of the three water layers and determining a transport of the particle comprises comparing the shear stress value in the water layer in which the particle is located to a shear stress threshold depending on the particle:

• • if the particle is located in the plume or subsurface layer, and,
    • if the shear stress value in the water layer in which the particle is located is higher than the suspension shear stress threshold, transporting the particle to an adjacent cell, said adjacent cell being determined based on the determined direction of the shear stress, or
    • if the shear stress induced on the particle is lower than the suspension shear stress threshold, transporting the particle towards the subsurface or bottom layer, respectively, within the same cell, or
  • if the particle is located in the bottom layer, and
    • if the shear stress induced on the particle is higher than a traction shear stress threshold, transporting the particle to an adjacent cell, said adjacent cell being determined based on the determined direction of the shear stress within the bottom layer, or
    • if the shear stress induced on the particle is lower than the traction shear stress threshold, depositing the particle.

In embodiments, the method comprises a preliminary step of determining a thickness of an Ekman layer extending from the surface of the immersed area, and defining the three water layers such that:

• • the bottom layer extends between the water bottom and a fixed distance thereof,
  • the plume layer extends between the water surface and a depth determined based on the Ekman layer's depth, and
  • the subsurface layer extends between the plume layer and the bottom layer.

In embodiments, the preliminary step further comprises defining parameters of wind occurring over the immersed area, and the method further comprises:

• • determining the direction and velocity of a wind-induced current occurring in the plume layer based on the wind parameters, and
  • determining the direction and velocity of a return current occurring in the subsurface layer and resulting from the wind-induced current occurring in the plume layer.

In embodiments, determining the direction and velocity of a wind-induced current occurring in the plume layer comprises determining the direction and velocity of an ocean surface current caused by an Ekman vortex.

In embodiments, determining the direction and velocity of a wind-induced current occurring in the plume layer comprises determining a direction and velocity of a wave-induced current.

In embodiments, determining the transport of a particle further comprises a preliminary modelling of aggregation of particles by flocculation, and the determining a transport of particles is implemented on the aggregated particles.

In embodiments, modelling aggregation of particles by flocculation comprises computing a proportion of aggregated particles based on a set of parameters including the computed shear stress and the particles size.

According to another aspect, a computer-implemented method of modelling sedimentary deposition within an immersed area is disclosed, comprising:

• • a setup step comprising defining a geological gridded model of the immersed area, setting a reference water level, at least one supply or production process of particles to be introduced within the model, and wind parameters, and
  • a step of simulating the evolution of the geological gridded model over a predetermined period of time T, comprising:
    • assigning a water depth to a plurality of cells,
    • determining a direction and velocity of at least one water current occurring within the immersed area,
    • determining, from the water current, a direction and intensity of a shear stress induced by the water current on a particle,
    • Introducing at least one particle in at least one cell of the geological gridded model,
    • determining a transport of the particle, from the determined direction and intensity of the shear stress induced on the particle, the granulometry and sediment type of the particle, and,
    • updating the geological gridded model of the area according to the transport of the particles.

In embodiments, the step of transporting the particles is repeated until all introduced particles are deposited or have exited the geological gridded model.

In embodiments, each supply or production process is chosen among one of the following groups:

• • clastic supply processes, comprising at least river mouth supply and mineral spring causing travertine deposition,
  • carbonates production process, and
  • each supply or production process is associated to a depth at which a particle is introduced.

According to another aspect, a computer program product is disclosed, comprising code instructions for implementing the methods according to the description above, when it is executed by a processor.

According to another aspect, a non-transitory computer readable storage medium is disclosed, having stored thereon a computer program comprising program instructions, the computer program being loadable into a processor and adapted to cause the processor to carry out, when the computer program is run by the processor, the methods according to the description above.

The method allows computing the transport induced by water currents on sedimentary particles, based on a shear stress induced by the water current on the particle. Relying on a shear stress rather than, for instance, a current velocity, allows differentiating the impact of a water current on sedimentary particles according to the depth of the water current. The result of the simulation is therefore rendered more realistic and can be more accurately compared with logs of existing reservoirs.

BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWINGS

The present disclosure is illustrated by way of example, and not by way of limitation, in the figures of the accompanying drawings, in which like reference numerals refer to similar elements and in which:

FIG. 1 is an exemplary representation of a geological gridded model of an area, showing a first time-layer and a cut view of the superposition of other time-layers above the first, where the grey level of the cells indicates the water depth of each cell,

FIG. 2 is a flow chart describing a possible embodiment of the method for modelling the water-current induced transport of particles,

FIG. 3a is a flow chart describing a possible embodiment of a method for modelling the evolution of a sedimentary area,

FIG. 3b is a flow chart describing another possible embodiment of the method for modelling the water-current induced transport of particles, taking into account occurrence of exceptional stormy events,

FIG. 4 is a representation of the decomposition of an immersed area into three water depths ranges,

FIG. 5 schematically represents the wind-induced currents and the evolution of their velocities according to water depth,

FIG. 6a represents a vertical profile of current velocity in the plume and subsurface layers,

FIG. 6b represents a vertical profile of shear stress in the plume and subsurface layers,

FIG. 7a represents the value of an avalanche angle as a function of grain size,

FIG. 7b represents the gravity-induced shear stress as a function of granulometry and topographic slope,

FIGS. 8a to 8d represent Shields diagrams showing whether a particle is deposited or transported according to its size and the shear stress applied to the particle, for four types of particles,

FIGS. 8e to 8h represent Shields diagrams showing whether a particle is remobilized according to its size and the shear stress applied to the particle, for four types of particles,

FIGS. 9a to 9d represent controlling factors for flocculation process,

FIG. 10 represents the distribution of particles before and after modelling a flocculation process, and

FIG. 11 is a possible embodiment for a device that enables the present disclosure.

DETAILED DESCRIPTION

The method described below models the evolution of sedimentary basins by simulating the deposition over time of clastic and/or carbonates particles that can be supplied or produced by diverse processes such as rivers, travertine sources, in-situ production of carbonates, remobilization of carbonates, etc. This method also takes into account the impact of water currents on the transport of clastic and carbonates particles. The sedimentary area to be modelled can be either a marine area, or a lacustrine area.

The method is a forward stratigraphic modelling method, modelling the evolution of an immersed sedimentary area through the deposition of successive layers of sediments. Accordingly, the sedimentary area is formed by a stack of layers, wherein each layer is defined by a two-dimensional surface representing the surface of the ground of the immersed area at the time of deposition of sediments, and by a thickness of sediments forming said layer. Each layer is gridded and comprises, like in the example shown in FIG. 1, a plurality of grid cells M_1,1, M_1,2, M_2,1, . . . , and more generally M_i,j, where the variables i and j indicate the positions of the cells within the surface. Each layer if further associated to a time t corresponding to a time of deposition of sediments. Accordingly, each cell of the modelled area can be described by three parameters (i, j, k) or (x,y,t) wherein the two first coordinates represent the location of the cell within the surface and k represents the time layer or t represents the time of formation of the surface, which is equivalent.

Typically, each cell represents an area having a side length of a few hundreds of meters, up to a few kilometers. The method starts with an initial topography corresponding to an initial surface representing the bottom of the immersed area, above which a column of water of defined height (i.e., water depth of each cell of the model) is defined, and comprises iterating a series of steps modelling the introduction into the model of clastic and/or carbonates particles during a predetermined period of time T, their transport induced by water currents during this period of time T, and the deposition of some of these particles to form an additional layer of sediments.

At the end of this period of time T, the topography of the model is updated in each cell according to the quantity of deposited particles. More specifically, a layer is generated, which thickness in each cell is determined based on a number of particles deposited at this cell. Such a layer is called a time layer since it corresponds to the passage of the predetermined period of time T. The topography of the geological gridded model of the area thus evolves with the accumulation of time layers.

The example shown in FIG. 1 represents the gridded model in which the main represented surface is a time-layer and a cross-section of the stack of superposed time-layers is also shown. The height of a cell represents the thickness of the deposited layer of sediments and hence the topography of the area. Moreover, the grey level of each cell represents its water depth at the considered time layer. In this example one can notice that the water-depth of the cells located on the right-hand side of the model is progressively reduced until being equal to 0 due to the successive deposition of sediments and hence the emersion of the corresponding land. Once the land is emerged the topography of the model no longer evolves since it is no longer subject to sedimentation.

The methods of modelling transport of particles and of modelling the evolution of a sedimentary area are computer-implemented methods. With reference to FIG. 11 is shown a possible embodiment for a device that enables the present disclosure.

In this embodiment, the device 10 comprises a computer, this computer comprising a memory 15 to store program instructions loadable into a circuit and adapted to cause circuit 14 to carry out the steps of the present disclosure when the program instructions are run by the circuit 14.

The memory 15 may also store data and useful information for carrying the steps of the present disclosure as described above.

The circuit 14 may be for instance:

• • a processor or a processing unit adapted to interpret instructions in a computer language, the processor or the processing unit may comprise, may be associated with or be attached to a memory comprising the instructions, or
  • the association of a processor/processing unit and a memory, the processor or the processing unit adapted to interpret instructions in a computer language, the memory comprising said instructions, or
  • an electronic card wherein the steps of the disclosure are described within silicon, or
  • a programmable electronic chip such as a FPGA chip (for «Field-Programmable Gate Array»).

This computer comprises an input interface 13 for the reception of several data used for the above method according to the disclosure, for instance the gridded model, some parameters of the topography of the modelled area, some parameters of the modelled currents, etc. This computer also comprises an output interface 16 for outputting the updated geological gridded model.

To ease the interaction with the computer, a screen 11 and a keyboard 12 or a tactile screen may be provided and connected to the computer circuit 14. The various components described above may be remotely connected to one another, i.e., the memory storing the data and/or the circuit implementing the method may be remotely located with reference to the user and accessible through any suitable network.

With references to FIGS. 3a and 3b the main steps of a method for modelling the evolution of a sedimentary area will now be disclosed. Among them, the main steps of a method for modelling water current-induced transport of particles are shown in FIG. 2.

Setup Step

The method for modelling the evolution of a sedimentary area comprises a first preliminary setup step 90.

The setup step comprises initializing a topography of the modelled area, i.e., receiving an initial surface having a plurality of cells wherein each cell corresponds to a position (x,y) and is assigned a parameter z0 which is the height of said bottom surface. The setup step also comprises defining the period of time T corresponding to one time layer, i.e., the period of time during which introduction of particles representing sediments, transport and deposition of these particles is simulated and, possibly, iterated. Preferably, this duration may be comprised between 1,000 and 10,000 years.

The setup step also comprises setting an initial reference water level z_r, as well as the evolution of the reference water level over the model between two successive periods of time T, i.e., two successive time layers (eustatism for marine areas) and the amount of subsidence of the ground's surface over the geological gridded model between two time layers. The amount of subsidence may vary over the model of the area, i.e., it may not be the same for all the cells of the models.

From the initial topography of the geological gridded model and the initial reference water level, a water depth WD in each cell is inferred and assigned to the respective cell. If the height z₀ of a cell is above the reference water level z_r, then the water depth is zero. It can be understood that as the method aims at modelling the evolution of a sedimentary area, at least some of the cells of the initial topography are below water level, i.e., z₀<z_r.

The setup step 90 also includes the user defining the number and types of particles supply or production processes, the number and types of water currents to be modelled, and the parameterizing of each particle supply or production process and of each current.

Regarding the particles supply or production processes, each particle introduced within the model may result from either a carbonate production model or a siliciclastic supply process.

Carbonates production processes comprise:

• • carbonates from underwater factories,
  • aerial travertine.

Siliciclastic supply processes comprise:

• • river mouth supply,
  • volcanoes supply.

The parameterization of production of carbonates may comprise defining a sedimentary element type, granulometry, and production rate of the carbonate production. The production rate may itself depend upon environmental factors which may also be set by the user, such as bathymetry, current energy, water temperature, chemistry, etc.

The parameterization of each siliciclastic supply process may comprise defining a location of the source of elements (location of the volcano, river mouth, etc.) sedimentary element type, granulometry, and rate of supply of the considered process, i.e., a volume or mass of supplied sedimentary elements per time unit.

During the subsequent modelling of the supply or production processes, the sedimentary elements are introduced in the model as particles where each particle represents a determined mass or volume of siliciclastic or carbonates sediments of a defined granulometric class.

Regarding the currents to be modelled, the user may choose to model at least one among the following water currents:

• • wind-induced current, including wind-induced wave current and oceanic surface current,
  • tidal current,
  • river-mouth induced current.

Oceanic surface current can only apply for marine areas, whereas the other currents, including tidal current, can apply for both marine and lake areas. If wind-induced currents are modelled, then the setup step also comprises setting up a celerity vector of the wind {right arrow over (u_wind)}, comprising setting a direction and an absolute value of the wind celerity.

If at least one river mouth induced current is to be modelled, the setup step 90 can comprise the user setting the volumetric flow of the river at the river mouth, the width and depth of the river mouth.

In embodiments, the modelling of evolution of a sedimentary area may also comprise simulating exceptional weather events, i.e., storms, floods, etc., that are prone to remobilizing sediments that have already been deposited. In these embodiments, the setup step 90 also comprises defining a duration of stormy conditions over the time period T, expressed for instance as a number of days of stormy conditions per year. If wind-induced currents are modelled, a celerity vector of the wind {right arrow over (u_wind,storm)} in stormy conditions may also be defined during this setup step. By contrast, in these embodiments, the celerity vector of the wind {right arrow over (u_wind)} introduced above relates to the parameters of the wind in fair weather conditions, i.e., during the rest of the days of the year.

The method then comprises a series of steps which are detailed below, and which are implemented to generate one time layer, representing the passage of the predetermined period of time T.

Computation Layer

The series of steps which forms one computation layer, and which is implemented at least once to generate a time layer representing the period of time T, or iterated a number N of times in the embodiments described above, is designated by reference 900 on FIGS. 3a and 3b, and will now be described.

An optional preliminary step 99 comprises the change of some parameters of the model by the user, if it is desired to represent an evolution of these parameters between one period of time represented by a time layer T and another. For instance, the parameters regarding the river mouth current that can be set at step 90 (volumetric flow, width and height at the river mouth) can be modified at step 99. Also, the parameterization of each supply or production process may be changed at optional step 99.

The eustatism and subsidence rate may also be amended between two time-layers during said preliminary step 99.

Step 100 comprises receiving the geological gridded model of the area, either by loading an initial version of the model, or by updating the model according to a previous iteration of the series of steps 900. The update comprises updating a height along z of each cell, which corresponds to the initial position z₀ of the cell along z added to the thickness of particles deposited at the cell. The height along z may also take into account local subsidence of the ground's surface.

The method then comprises a step 200 of computing, from the topography updated at step 100, topographic slopes of the ground surface formed by the model and inferring, from the topography and reference water level z_r, the water depth WD in each cell.

The method then comprises a step 300 of modelling marine water currents occurring over the immersed area represented by the gridded model. This step is performed by determining, for a plurality of cells of the gridded model for which WD>0, and preferably each cell for which WD>0, a direction and velocity of each water current to be modelled. The shoreline SL is defined by cells for which WD=0 and which are adjacent cells for which WD>0.

For the computation of direction and velocity of each water current, and with reference to FIG. 4, the extent of water extending over the ground surface is decomposed into three water layers or respective depths, comprising a bottom layer, extending at water bottom, a plume layer, extending at water surface, and a subsurface layer extending between the plume layer and the bottom layer.

More specifically, the thickness of the bottom layer is constant and corresponds to the boundary layer for the subsurface flow. It is preferably comprised between a few mm and a few cm. It may be set by the user or defined per default. According to a non-limiting example, the thickness of the bottom layer may be set equal to 0.01 meter. The bottom layer thus extends between the ground surface and a fixed distance thereof.

The thickness of the plume layer is determined based on the thickness of Ekman layer.

Ekman layer is the layer in a mass of fluid on Earth where there is a force balance between pressure gradient force, Coriolis force and turbulent drag. The thickness of the Ekman layer depends upon the latitude and is expressed as follows:

$H_{\mathrm{Ekman}}=\{\begin{array}{c}\pi \sqrt{2A_z/❘f❘}\mathrm{if}\varphi \ge 7.7°\\ 100\mathrm{if}\varnothing \le 7.7°\end{array}$

where A_z=0.01 m²/s, f=Ω·sin(Φ) and Q=7.3·10⁵ rad/s.

The thickness of the plume layer, i.e., the limit between the plume and subsurface layer, is set to be at most equal to H_Ekman computed above if the water depth in a given cell is at least 2H_Ekman. Otherwise, the thickness of the plume layer is equal to WD/2 as summarized below:

$z_{\mathrm{plume}/\mathrm{subsurf}}=\{\begin{array}{cc}{z}_{\mathrm{surface}}-H_{\mathrm{Ekman}}& \mathrm{if}\mathrm{WD}\ge 2H_{\mathrm{Ekman}}\\ z_{\mathrm{surface}}-\frac{\mathrm{WD}}{2}& \mathrm{if}\mathrm{WD}\le 2H_{\mathrm{Ekman}}\end{array}$

An example of decomposition into three water layers depending on the water depth of the ground surface and the thickness of the Ekman layer is shown in FIG. 4.

Computing 300 the direction and velocity of currents occurring within the immersed area thus comprises determining 310 the thickness of the Ekman layer, decomposing the water into three water layers 320, and then determining 330 a direction and velocity of currents occurring within the plume layer, and within the subsurface layer. In the bottom layer, the velocity of currents is not computed since it is considered as the boundary layer of the subsurface current.

As indicated above, one or more currents may be simulated among: wind-induced currents, which comprise ocean surface current and wave-induced current, river-mouth current and tidal current. Accordingly, a celerity vector may be determined for each current in the plume and subsurface layer.

The celerity vector of the wind-induced currents in the plume layer {right arrow over (U_WIC,plume)} is the sum of the celerity vectors of the ocean surface current {right arrow over (U_OSC,plume)} and the wave-induced current {right arrow over (U_waves,plume)}:

$\overrightarrow{U_{\mathrm{WIC},\mathrm{plume}}}=\overrightarrow{U_{\mathrm{OSC},\mathrm{plume}}}+\overrightarrow{U_{\mathrm{waves},\mathrm{plume}}}$

Ocean surface current is the current induced in the water surface by the wind, inducing a water displacement within the Ekman layer.

In the plume layer, the direction of displacement of water induced by the ocean surface current (hereinafter denoted “OSC” in indexes) is perpendicular to the direction of the wind, towards the right (i.e., clockwise direction) in the Northern hemisphere and towards the left (i.e., counter-clockwise direction) in the Southern hemisphere:

$\overrightarrow{U_{\mathrm{OSC},\mathrm{plume}}}=\{\begin{array}{c}\begin{array}{c}\overrightarrow{U_{{\mathrm{OSC}}_0}}=\pm \frac{{\rho }_{\mathrm{air}}{C}_f{❘u_{\mathrm{wind}}❘}^2}{{\rho }_w.f.H_{\mathrm{Ekman}}}.\left(-{\mathrm{Dir}}_{{\mathrm{wind}}_y},{\mathrm{Dir}}_{{\mathrm{wind}}_x}\right)\mathrm{if}\\ \mathrm{WD}\ge 2H_{\mathrm{Ekman}}\end{array}\\ \overrightarrow{U_{{\mathrm{OSC}}_0}}.\frac{\mathrm{WD}}{2H_{\mathrm{Ekman}}}\mathrm{if}\mathrm{WD}\le 2H_{\mathrm{Ekman}}\end{array}$

where Dir_windx and Dir_windy are the coordinates in x and y of the wind speed vector, ρ_w and ρ_air are the volumetric masses of water and air, respectively, in kg/m³, C_f=0.0015 is the friction coefficient of water surface. The norm of the vector U_OSC,plume is the velocity value of the OSC current in the plume layer in a considered cell. It is an average value over the thickness of the plume layer.

In the plume layer, the wave-induced current is parallel to the direction of the wind and only occurs from the surface down to a determined depth WD_base called wavebase water depth, so by definition this wave-induced current only occurs in cells in which the water depth WD is below WD_base. WD_Base may be computed as follows:

${\mathrm{WD}}_{\mathrm{Base}}=\frac{{\lambda }_{\mathrm{wave}}}{2}$

where λ_wave is the wave wavelength induced by the wind and may be inferred from the wind speed by use of a wave model such as the Airy linear model.

A wave breaking water depth is also defined as:

${\mathrm{WD}}_{\mathrm{wavebreaking}}=\frac{{\lambda }_{\mathrm{wave}}}{20}$

When stormy events are simulated these parameters have to be computed for both stormy and fair weather conditions.

The velocity of the wave-induced current is computed from the height of the waves and the celerity of the waves.

The height of the waves is computed as follows:

$H_{\mathrm{waves}}=\{\begin{array}{c}{H}_{{\mathrm{waves}}_0}=\frac{2C_{{\mathrm{waves}}_0}^2}{9g}\mathrm{if}\mathrm{WD}={\mathrm{WD}}_{\mathrm{base}}\\ H_{{\mathrm{waves}}_0}.\sqrt{\frac{C_{{\mathrm{waves}}_0}}{C_{\mathrm{waves}}}}\mathrm{if}{\mathrm{WD}}_{\mathrm{wavebreaking}}\le \mathrm{WD}\le {\mathrm{WD}}_{\mathrm{base}}\\ \begin{array}{c}\mathrm{linear}\mathrm{interpolation}\mathrm{between}{H}_{{\mathrm{waves}}_0}.\sqrt{\frac{C_{{\mathrm{waves}}_0}}{C_{\mathrm{waves}}}}\mathrm{and}\\ 0\mathrm{if}\mathrm{WD}\le {\mathrm{WD}}_{\mathrm{wavebreaking}}\end{array}\end{array}$

The celerity of the waves is computed as follows:

$C_{\mathrm{waves}}10=\{\begin{array}{c}{C}_{{\mathrm{waves}}_0}=U_{\mathrm{wind}}\mathrm{if}\mathrm{WD}={\mathrm{WD}}_{\mathrm{base}}\\ \sqrt{g.\mathrm{WD}}\mathrm{if}\mathrm{WD}\le {\mathrm{WD}}_{\mathrm{wavebreaking}}\\ \begin{array}{c}\mathrm{Linear}\mathrm{interpolation}\mathrm{between}{C}_{{\mathrm{waves}}_0}\mathrm{and}\\ \sqrt{g.{\mathrm{WD}}_{\mathrm{wavebreaking}}}\mathrm{if}{\mathrm{WD}}_{\mathrm{wavebreaking}}\le \mathrm{WD}\le {\mathrm{WD}}_{\mathrm{base}}\end{array}\end{array}$

The velocity vector of the wave induced current is thus null if WD>WD_base and is equal, when WD≤WD_base to:

$\overrightarrow{U_{\mathrm{waves},\mathrm{plume}}}=\frac{{\mathrm{gH}}_{\mathrm{waves}}^2}{8.\mathrm{WD}.C_{\mathrm{waves}}}\overrightarrow{{\mathrm{Dir}}_{\mathrm{wind}}}$

The norm of the vector U_waves,plume is the velocity of the wave-induced current in the plume layer in a considered cell. It is an average velocity over the thickness of the plume layer. Additionally, when the water depth WD is below twice the thickness of the bottom water layer, the velocity of the wave-induced current is set to zero to prevent numerical divergence.

FIG. 5a provides a summary of the ocean surface current velocity and wave induced current velocity in the plume layer.

In embodiments, the modelling of wind-induced currents may also comprise modelling a fetch distance. Indeed, the waves are formed only when the wind blows on wide enough areas. Thus when the wind blows from the sea towards the shore line, it is considered that the wave-induced current and oceanic surface current appear in all the cells until the shoreline according to the above equations. When by contrast the wind blows from the shoreline towards the sea, a fetch distance is computed from the shoreline, which value is defined by:

Fetch=0.20·U_wind^0.57

If the water depth of the cells is lower than WD_base from the shoreline to the fetch distance, the celerity of the waves linearly increases from 0 on the shoreline until Uwaves computed according to the equation above at the fetch distance. If the water depth is higher than WD_base then the wave celerity is null.

Regarding the ocean surface current, the latter increases linearly from 0 at the shoreline until reaching the velocity U_OSC computed according to the equation above at the fetch distance.

Regarding wind-induced currents, the direction and velocity of these currents in the subsurface layer can be inferred at substep 340 from the direction and velocity respectively of the currents in the plume layer.

In the subsurface layer, with reference to FIG. 5b, the wind-induced currents are return currents (denoted {right arrow over (U_{return,subsurface})}) from the wind-induced currents occurring in the plume layer.

Accordingly, the velocity vectors of said return currents extend parallel to the topographic slope in the considered cell, i.e., the vector of maximum slope of the ground surface in the cell, and in a direction opposed to the direction of the component of the respective velocity vector of the wind induced currents within the plume layer that is parallel to the slope. It is assumed that the component orthogonal to the slope of the wind induced currents within the plume layer—i.e., longshore component—keeps flowing away from the considered cell.

The subsurface velocity vector is thus inferred from the component of the velocity of the wind-induced currents in the plume layer which is parallel to topographic slope as follows,

$\overrightarrow{U_{\mathrm{return},\mathrm{subsurface}}}=-\left(\overrightarrow{U_{\mathrm{WIC}}}.\overrightarrow{\mathrm{Slope}}\right).\overrightarrow{\mathrm{Slope}}.\frac{H_{\mathrm{plume}}}{H_{\mathrm{subsurface}}}$

{right arrow over (Slope)} is the unitary vector of maximum descending slope in the considered cell.

If, in addition to the wind-induced currents, the modelled currents also comprise a river-mouth induced current and/or a tidal current, then:

• • the directions and velocities of the river-mouth induced current in the plume and in the subsurface layers may also be determined by implementing the method disclosed in WO 2020/229865.

In particular, the velocity of a river-mouth induced current depends on flow type of the river and the position of the considered cell with respect to the width of the river jet. The river jet type may be either homopycnal, hypopycnal and hyperpycnal and may be determined based on a comparison between the density of the river jet (water and sedimentary discharge) and the density of the water in which the river flows. The velocity vector of river-mouth induced current extends perpendicular to the direction along which extends the width of the river mouth, and extending away from the river mouth, and the velocity of the current U_river is denoted as u in WO2020/229865.

• • the direction and velocity of the tidal current may also be determined by implementing the method disclosed in WO 2020/229866.

In particular, as developed in this document, modelling the tidal current implies decomposing each time layer into an even number of 2k computation layers corresponding to subperiods of time of duration T/2k where half of the subperiods correspond to rising tide and the other half corresponds to the falling tide. The direction of the tidal current is towards the shoreline during rising tide and away from the shoreline during falling tide. The velocity of the tidal current U_tidal in the plume layer increases from 0 at the high tide shoreline to a maximum value (which may be user-defined) at a distance from the shoreline such that WD=z_r in the cell, and then decreases again until reaching 0 in cells where WD=WD_TD where WD_TD is the water depth of influence of the tidal current. The velocity of the tidal current in the subsurface layer is computed from the velocity of the tidal current in the plume layer of the same cell by applying a decrease factor which is function of water depth.

The method then comprises determining 400 a direction and intensity of shear stress induced by the modelled water currents in each water layer. The intensity of a shear stress in a water layer is a depth-averaged intensity over the height of said water layer. It is thereafter denoted:

• • {right arrow over (τ_plume)} a shear stress induced by simulated currents on particles in the plume layer,
  • {right arrow over (τ_subsurface)} a shear stress induced by simulated currents on particles in the subsurface layer, and
  • {right arrow over (τ_bottom)} a shear stress in the bottom layer, which is induced as explained below by the simulated currents and also the slope of the bottom. The part of shear stress in the bottom layer only induced by water currents is denoted {right arrow over (τ_bottom,cur)}.

The shear stress in a water layer is the sum of the shear stresses induced in that water layer by the simulated currents. For instance, when the simulated currents include wind-induced currents, river-mouth current and tidal current, the shear stress in the plume layer is computed as follows:

$\overrightarrow{{\tau }_{\mathrm{plume}}}=\overrightarrow{{\tau }_{\mathrm{WIC},\mathrm{plume}}}+\overrightarrow{{\tau }_{\mathrm{river},\mathrm{plume}}}+\overrightarrow{{\tau }_{\mathrm{tidal},\mathrm{plume}}}\mathrm{where}\overrightarrow{{\tau }_{\mathrm{tic},\mathrm{plume}}}\left(z\right)=\overrightarrow{{\tau }_{\mathrm{waves},\mathrm{plume}}}\left(z\right)+\overrightarrow{{\tau }_{\mathrm{OSC},\mathrm{plume}}}\left(z\right)$

The depth-average shear stress induced by a water current in a water layer is determined from a shear stress profile according to depth, which is itself determined from a velocity profile of the water currents as a function of water depth in said layer.

Regarding the wind-induced currents, i.e., the wave-induced current and the ocean surface current, each velocity profile in the plume layer is a power law profile, wherein the velocity is maximum at the surface and becomes null at the interface between the plume and the subsurface layers, such that:

$u_{\mathrm{waves}}\left(z\right)={U_{\max ,\mathrm{waves}}\left(\frac{z}{H_p}\right)}^m{u}_{\mathrm{OSC}}\left(z\right)={U_{\max ,\mathrm{OSC}}\left(\frac{z}{H_p}\right)}^{1/m}$

With m is the power of the velocity profile, here m=4, and H_p is the thickness of the plume layer. The maximum values U_max,OSC and U_max,waves at the surface may be determined by computing the integral from the surface to the interface between plume and subsurface layer of each velocity profile as follows:

$U_{\max ,\mathrm{OSC}}=U_{\mathrm{OSC}}\frac{1+m}{m}{U}_{\max ,\mathrm{waves}}=U_{\mathrm{waves}}\left(1+m\right)$

where U_OSC and U_waves are the mean velocities, over the thickness of the plume layer, of each current, computed in step 330 above.

The shear stress profile according to depth within the plume layer can be computed based on the velocity profile by computing a derivative with respect to depth of the velocity profile of the currents:

$\tau \left(z\right)=k\frac{\partial u\left(z\right)}{\partial z},$

where k is a calibration constant, for instance equal to 0.1.

${\tau }_{\mathrm{waves},\mathrm{plume}}\left(z\right)={{\tau }_{0,\mathrm{waves}}\left(\frac{z}{H_p}\right)}^{m-1}{\tau }_{\mathrm{OSC},\mathrm{plume}}\left(z\right)={{\tau }_{0,\mathrm{OSC}}\left(\frac{z}{H_p}\right)}^{1/m-1}\mathrm{where}{\tau }_{0,\mathrm{waves}}=U_{\mathrm{waves}}.\frac{k.m.\left(1+m\right)}{H_p}\mathrm{and}{\tau }_{0,\mathrm{OSC}}=U_{\mathrm{OSC}}.\frac{k.\left(1+m\right)}{H_p.m^2}.$

Once the shear stress profile is obtained, computing its average along the vertical axis enables obtaining a mean value of the shear stress induced by the wind-induced currents within the plume layer. Further, the directions of the shear stresses induced by the wave-current and the ocean surface currents, respectively, is the same as the direction of the respective velocity vectors

$\overrightarrow{{\tau }_{\mathrm{wic},\mathrm{plume}}}=\overrightarrow{{\tau }_{0,\mathrm{OSC}}}.m+\overrightarrow{{\tau }_{0,\mathrm{waves}}}.\frac{1}{m}$

Regarding computation of a shear stress value incurred by the wind-induced currents within the subsurface and bottom layers, as these currents extend in an opposite direction in the subsurface layer and in the plume layer, it is considered that the velocity of these currents is null at the interface plume/subsurface layer (z=H_s=H_sub+H_bot where H_bot is the thickness of the bottom layer and H_sub is the thickness of the subsurface layer) and also null on the ground, and the velocity reaches a maximum between those two null points. The velocity profile in the subsurface and bottom layers u_sub(z) can thus be written as a sum of velocity profiles occurring within the subsurface and bottom layers, where the velocity profiles are power law profiles of opposite signs with respectively m and 1/m exponents: u_sub(z)−u₁(z)+u₂(z)

$u_1\left(z\right)={U_{\max ,1}\left(\frac{z}{H_s}\right)}^m{u}_2\left(z\right)={U_{\max ,2}\left(\frac{z}{H_s}\right)}^{1/m}$

where u₁(z) designates is a first velocity profile of this decomposition, U_max,1 is the corresponding maximum velocity, u₂(z) designates the second velocity profile of this decomposition, and U_max,2 is the corresponding maximum velocity.

Computing the integral between the plume/subsurface boundary and the ground of the velocity profile, and considering the null boundary conditions, allows obtaining the two constant parameters U_max,1 and U_max,2:

$U_{\max ,2}=-U_{\max ,1}=U_{\mathrm{wic},\mathrm{sub}}\frac{1+m}{m-1}$

where U_wic,sub is the mean velocity of return currents in the subsurface layer computed above.

The shear stress profile according to depth within the subsurface and bottom layers of the return currents can be computed by computing the derivative of the current velocity profile in the subsurface layer:

$\overrightarrow{{\tau }_{\mathrm{sub}}}\left(z\right)=\overrightarrow{{\tau }_1}\left(z\right)+\overrightarrow{{\tau }_2}\left(z\right)=\to \overrightarrow{{\tau }_{0,1}}·{\left(\frac{z}{H_s}\right)}^{m-1}+\overrightarrow{{\tau }_{0,2}}·{\left(\frac{z}{H_s}\right)}^{1/m-1}$ $\mathrm{With}\overrightarrow{{\tau }_{0,1}}=\overrightarrow{U_{sub}}·\frac{k·m·\left(1+m\right)}{H_s\left(m-1\right)}\mathrm{and}\overrightarrow{{\tau }_{0,2}}=-\overrightarrow{U_{sub}}·\frac{k·\left(1+m\right)}{H_s\left(m-1\right)·m}.$

Once the shear stress profile is obtained, the mean shear stress induced by return currents within the subsurface layer, i.e., from the interface between the bottom and subsurface layers (z=H_bot) and the interface between the subsurface and plume layers (z=H_s), can be computed as:

$\overrightarrow{{\tau }_{\mathrm{wic},\mathrm{subsurface}}}=\{\begin{array}{cc}\frac{1}{1-z_b}\left(\frac{\overrightarrow{{\tau }_{0,1}}}{m}\left(1-2a^m+z_b^m\right)+m·\overrightarrow{{\tau }_{0,2}}\left(1-2a^{\frac{1}{m}}+z_b^{\frac{1}{m}}\right)\right)& \mathrm{if}a\ge z_b\\ \frac{1}{1-z_b}\left(\frac{\overrightarrow{{\tau }_{0,1}}}{m}\left(1-z_b^m\right)+m·\overrightarrow{{\tau }_{0,2}}\left(1-z_b^{\frac{1}{m}}\right)\right)& \mathrm{if}a\le z_b\end{array}$ $\mathrm{With}{z}_b=\frac{H_{bot}}{H_s}\mathrm{and}a={\left(\frac{1}{m}\right)}^{\frac{m}{m^2-1}}.$

The velocity profiles of the wind-induced currents within the plume and subsurface layers are shown in FIG. 6a, and the shear stress profiles induced by said currents within the same water layers are shown in FIG. 6b.

As the bottom layer corresponds to the boundary limit of the subsurface layer's flow, a mean shear stress value in said layer is computed by the mean of the shear stress profile between the ground (z=0) and the top of the bottom layer, i.e., the interface between the bottom and subsurface layers (z=H_bot), and is equal to (considering only return currents of wind-induced currents):

$\overrightarrow{{\tau }_{\mathrm{wic},\mathrm{bottom}}}=\{\begin{array}{cc}\frac{1}{z_b}\left(\frac{\overrightarrow{{\tau }_{0,1}}}{m}{z}_b^m+m·\overrightarrow{{\tau }_{0,2}}{z}_b^{\frac{1}{m}}\right)& \mathrm{if}a\ge z_b\\ \frac{1}{z_b}\left(\frac{\overrightarrow{{\tau }_{0,1}}}{m}\left(z_b^m-z_a^m\right)+m·\overrightarrow{{\tau }_{0,2}}\left(z_b^{\frac{1}{m}}-2a^{1/m}\right)\right)& \mathrm{if}a\le z_b\end{array}$

Regarding the tidal current and river-mouth current, their velocity profile according to water depth is also assumed to be a power function with a 1/m exponent, but said velocity profile extends over all the water depth from surface to bottom:

$u_{\mathrm{river}}={U_{\max ,river}\left(\frac{z}{\mathrm{WD}}\right)}^{1/m}$ $u_{\mathrm{tidal}}\left(z\right)={U_{\max ,\mathrm{tidal}}\left(\frac{z}{\mathrm{WD}}\right)}^{1/m}$

The maximum values U_max,river and U_max,tidal at the surface may be determined by computing the integral from the surface to the interface between plume and subsurface layer of each velocity profile:

$U_{\max ,\mathrm{river}}=U_{\mathrm{river}}\frac{1+m}{m}$ $U_{\max ,\mathrm{tidal}}=U_{\mathrm{tidal}}\frac{1+m}{m}$

U_river is equal to a defined in WO 2020/229865. U_tidal corresponds to the value of velocity E_tc,plume of the tidal current defined in WO 2020/229866.

The shear stress profile is then a derivative of the velocity profile according to z:

$\overrightarrow{{\tau }_{\mathrm{river}}}\left(z\right)={\overrightarrow{{\tau }_{0,r\mathrm{iver}}}\left(\frac{Z}{\mathrm{WD}}\right)}^{1/m-1}$ $\overrightarrow{{\tau }_{t\mathrm{idal}}}\left(z\right)={\overrightarrow{{\tau }_{0,t\mathrm{idal}}}\left(\frac{Z}{\mathrm{WD}}\right)}^{1/m-1}$ $\mathrm{where}\overrightarrow{{\tau }_{0,ri\mathrm{ver}}}=\overrightarrow{U_{r\mathrm{iver}}}·\frac{k\left(1+m\right)}{\mathrm{WD}·m^2}\mathrm{and}\overrightarrow{{\tau }_{0,ti\mathrm{dal}}}=\overrightarrow{U_{t\mathrm{idal}}}·\frac{k\left(1+m\right)}{\mathrm{WD}·m^2}$

The shear stress vector in each water layer can then be inferred by this profile by integrating the shear stress profile between the limits of each current, and the direction of the vector is provided by the direction of the velocity vector of the corresponding current.

Finally, the total shear stress vectors in the plume and in the subsurface layers are the sums of the shear stress vectors resulting from each simulated current in the corresponding layer.

Additionally, a specificity of the bottom layer is that the sediments are not only submitted to a shear stress induced by the flowing currents, but also to gravity, leading to take into account a gravity component of the shear stress within the bottom layer. This component depends on the topographic slope of the ground surface but also on proprieties of the sediments. Accordingly:

• • if the slope is null, the shear stress induced by the gravity is null,
  • the value of the shear stress induced by gravity on the particle is equal to a deposition shear stress threshold when the current velocity in the subsurface layer is null and the topographic slope of the ground surface is superior or equal to an avalanche angle determined for the particle, and
  • the direction of the shear stress induced by gravity on the particle is parallel to the direction of the downward topographic slope of the water bottom.

The gravity-induced shear stress applied to a particle of granulometry gr is computed as follows:

$\overrightarrow{{\tau }_{\mathrm{slope}.\mathrm{gr}}}=\frac{1}{e-1}\left(\exp \left(\frac{\alpha }{{\alpha }_{\mathrm{cr}}\left(\mathrm{gr}\right)}\right)-1\right)·{\tau }_{\mathrm{Cr}}·\overrightarrow{\mathrm{slope}}$

where α is the topographic slope, αc_r is the avalanche angle, i.e., the critical slope of the particles, depending on their granulometry as shown in FIG. 7a, τ_Cr is the critical shear stress of deposition of the particles, and {right arrow over (slope)} is the direction of the descending slope. FIG. 7b shows the gravity-induced shear stress as a function of granulometry.

The total shear stress within the bottom layer applied to particle of granulometry gr is therefore equal to:

$\overrightarrow{{\tau }_{bottom}}=\overrightarrow{{\tau }_{bott\mathrm{om},\mathrm{Cur}}}+\overrightarrow{{\tau }_{\mathrm{slope},\mathrm{gr}}}$

Back to FIGS. 3a and 3b, once the shear stress has been computed, the method further comprises a step 500 of introducing at least one particle of sediments in at least one cell of the geological gridded model. The number of particles introduced at step 500, the type of elements and their granulometric class, depend on the supply or production model defined at setup step 90.

The location, within the model, where the particles are introduced depends on the source of the particles. For instance, regarding clastic supply processes such as river mouth or volcano, the location of the source of particle is set during the setup step. Regarding in-situ production of carbonates, said production may occur in any cell where the conditions enabling the production are met. For instance, if a production function is associated with conditions on water depth, the production may occur only in the cells of the model in which the water depth satisfies the conditions.

The depth at which the particles are introduced also depend from the process they originate from and the granulometric class of the particle. Accordingly, particles originating from:

• • mineral sources causing deposit of travertine, or
  • in-situ production of carbonates, are introduced at the water bottom, immediately deposited, and may be remobilized when stormy events are simulated, as explained in greater details below. On the other hand, particles originating from a river mouth can be introduced in various water layers as detailed in WO2020/229865.

The method then comprises a step 600 of determining a transport of the particles of siliciclastic sediments introduced in the geological gridded model, based on the water-currents induced shear stress. The transport may comprise displacing a particle from a cell to a neighboring cell and/or depositing the particle.

FIGS. 8a to 8d represent Shields diagrams for deposition of particles, expressing the shear stress values for transport and deposition of sediments according to their grain diameter. FIG. 8a relates to porous flat carbonates, FIG. 8b relates to porous spherical carbonates, FIG. 8c relates to non porous spherical carbonates, and FIG. 8d relates to clast. The transport of a particle may be either a transport by traction or a transport by suspension. In FIGS. 8a to 8d, the area under the solid line corresponds to a zone where the particles are not transported, the area between the solid line and the dotted line corresponds to a zone where the particles are transported by traction, and the area above the dotted line corresponds to a zone where the particles are transported by suspension. The solid line thus corresponds to a traction shear stress threshold, and the dotted line corresponds to a suspension shear stress threshold. The determination of the transport mode of a particle depends on the shear stress induced on the particle at the water level where the particle is located and on the shear stress thresholds which depend on the particle, and which are represented by the lines of the Shields diagram.

Accordingly, the step of transporting and/or depositing a particle is performed as follows:

• • If the particle is located in the plume or subsurface layer, and
    • if the shear stress value in the water layer in which the particle is located is higher than the suspension shear stress threshold, the particle is transported to an adjacent cell, said adjacent cell being determined based on the determined direction shear stress, and the particle remains within the same water layer, or
    • if the shear stress induced on the particle is lower than the suspension shear stress threshold, the particle is transported downwards, i.e., to the subsurface layer if the particle is originally in the plume layer, or in the bottom layer if the particle is originally in the subsurface layer, while remaining within the same cell; or
  • if the particle is located in the bottom layer, and
    • if the shear stress induced on the particle is higher than a traction shear stress threshold, the particle is transported to an adjacent cell, said adjacent cell being determined based on the determined direction of the shear stress within the bottom layer, and remains within the bottom layer, or
    • if the shear stress induced on the particle is lower than the traction shear stress threshold, the particle is deposited on the ground surface.

In embodiments, the transport step 600 may comprise a preliminary substep 601 of modelling aggregation of particles by flocculation. Flocculation is a process by which fine-grain particles are electrically charged, which lead them to being gathered together due to Van des Waals forces. If flocculation is modelled, then preferably every implementation of the transport step includes said preliminary substep of modelling flocculation.

Accordingly, a fraction of flocculated particles is computed based on a plurality of controlling factors including at least the particle size and the current-induced shear stress occurring within the considered layer of water. The controlling factors preferably also comprise salinity and the concentration of suspended sediments, the latter being computed in each cell as V_particle/V_water where V_particle is the total volume of particles present in the cell (resulting from introduction and/or transport of particles into the cell) and V_water is the volume of water in the cell, which is computed from the side length of the cell and the water depth in the cell.

In embodiments, the fraction of flocculated particles may be computed based on the following function:

$F_{\mathrm{floc}}=T·\left(F_{\mathrm{floc},0}+\left(1-F_{\mathrm{floc},0}\right)·\prod _{\mathrm{controlling}\mathrm{factors}}{F}_i\right)$

Where F_floc,O is a background flocculation fraction, i.e., a proportion of particles which flocculate no matter the values of the controlling factors. It is a constant value which may be user-defined, and for instance set to 0. Fi is a controlling factor function expressing a flocculation fraction in % according to a value of the corresponding controlling factor.

With reference to FIGS. 9a to 9d are shown exemplary controlling factors functions. As shown in FIG. 9a, the fraction of flocculated particles decreases with the size of the particles, but increases with salinity (FIG. 9b) and with the concentration of suspended sediments (FIG. 9c). The fraction of flocculated particles also increases with the current induced layer shear stress until the latter reaches a threshold value (0.36 in FIG. 9d), and then decreases with the shear stress value. The shear stress value is here used as an indicator of the turbulence of the water: a turbulence increase brings together the particles and fosters contact between them; but when turbulence exceeds a threshold the particles no longer aggregate one another.

Flocculation is performed by aggregating the particles two by two starting from the finer ones, until the amount of aggregated particles reaches the fraction F_floc. For a considered particle, the particle that is aggregated to it is randomly cast according to the granulometric distribution of the particles. The result is a flocculated particle which size is the sum of the sizes of the two aggregated particles.

On FIG. 10 is shown an example of initial distribution of introduced particles at step 500 and the result of the flocculation process on said distribution.

Step 600 of transporting and/or depositing particles is iterated until all particles are deposited or have exiting the model.

When the modelling method includes simulation of stormy events, the method may further comprise a step 700 of modelling remobilization of a fraction of the particles deposited at step 600, following occurrence of stormy events. During implementation of step 700, only a fraction of the particles deposited during the same computation layer, corresponding to the same time-period, are remobilized. The particles deposited during a previous computation layer corresponding to a previous time-period are not remobilized.

Modelling remobilization 700 of a fraction of deposited particles comprises determining 710 a shear stress induced on the deposited particles during a stormy event and determining a fraction of the deposited particles which is remobilized, as well as the water current 720 in which the particles are remobilized particles.

Determining the shear stress induced on the deposited particles during a stormy event 710 involves determining a velocity and a direction of at least one wind-induced current occurring within the modelled area during a stormy event. This step can be performed according to the description of step 300 above, except that the celerity vector of the wind is the vector {right arrow over (u_w,storm)} assigned to stormy event, which is different from the celerity vector of the wind {right arrow over (u_w)} in fair weather conditions, i.e., the speed is greater and the direction may vary. When other currents are modelled during step 400, the same currents may also be modelled for stormy events, and a resulting current may be computed from the sum of modelled currents.

A shear stress value and direction imparted by the currents on the deposited particles, i.e., a shear stress value imparted by the currents in the bottom layer of the cells, is then computed in the same manner as step 400 disclosed above.

The fraction of remobilized particles of a given grain size gr is then computed as follows:

${\mathrm{Frac}}_{\mathrm{remobilized}}^{\mathrm{gr}}=\{\begin{array}{cc}0& \mathrm{if}{\tau }_{\mathrm{Bottom}}\le {\tau }_{\mathrm{threshold},\min }^{\mathrm{gr}}\\ t_f·\frac{{\tau }_{\mathrm{Bottom}}-{\tau }_{\mathrm{threshold},\min }^{\mathrm{gr}}}{{\tau }_{\mathrm{threshold},\max }^{\mathrm{gr}}-{\tau }_{\mathrm{threshold},\min }^{\mathrm{gr}}}& \mathrm{if}{\tau }_{\mathrm{threshold},\min }^{\mathrm{gr}}\le {\tau }_{\mathrm{Bottom}}\le {\tau }_{\mathrm{threshold},\max }^{\mathrm{gr}}\\ t_f& \mathrm{if}{\tau }_{\mathrm{Bottom}}\ge {\tau }_{\mathrm{threshold},\max }^{\mathrm{gr}}\end{array}$

where τ_{threshold,min}^gr and τ_{threshold,max}^gr are shear stress thresholds of the deposited particles which may be defined by the user. t_f is a time factor defined as:

$t_f=\{\begin{array}{cc}\frac{t_{\mathrm{storm}}}{t_{\max }}& \mathrm{if}{t}_{\mathrm{storm}}\le t_{\max }\\ 1& \mathrm{if}{t}_{\mathrm{storm}}>t_{\max }\end{array}$

And t_max is a parameter from which the time factor t_f is equal to 1. t_max can for instance be set at 10 days/year.

If the shear stress induced by the currents in the bottom layer is lower than the minimum shear stress threshold value τ_{threshold,min}^gr, then no particles are remobilized.

When particles are indeed remobilized (Frac_remobilized≠0), the water layer in which the particles are remobilized is determined 720 based on the shear stress value induced by the currents in at least the bottom layer—possibly in both the bottom and the subsurface layers, and possibly in each layer—and at least one characteristic shear stress value of the considered particle:

• • the fraction of remobilized particle is integrally remobilized in the bottom layer if the shear stress value induced by the currents in the bottom layer is greater than the critical motion shear stress value of the particles, but if the shear stress value induced by the currents in the subsurface layer is lower than a critical suspension shear stress value of the particles:

if τ_BOTTOM>τ_Cr,motion and τ_SUBSURFACE≤τ_{Cr,suspension}

then 100% of Frac_remob is remobilized in bottom layer

• • if the shear stress value induced by the currents in the subsurface layer is greater than the critical suspension shear stress value of the particles, and the shear stress value induced by the currents in the plume layer is lower than said critical suspension shear stress value, a part of the remobilized fraction of particles is remobilized in the subsurface layer and the rest is remobilized in the bottom layer:

if τ_BOTTOM>τ_Cr,motion and τ_SUBSURFACE>τ_{Cr,suspension} and τ_PLUME≤τ_{Cr,suspension}

then a suspended fraction F_susp is remobilized in the subsurface layer

• • if the shear stress values induced by the currents in both the subsurface and the plume layer are greater than the critical suspension shear stress value of the particles, then part of the remobilized suspended fraction of particles is remobilized in the plume layer and the rest is remobilized in the subsurface layer:

if τ_BOTTOM>τ_Cr,motion and τ_SUBSURFACE>τ_{Cr,suspension} and τ_PLUME>τ_{Cr,suspension}

then a fraction F_plume of the suspended fraction F_susp is remobilized in the plume layer and 1−F_plume is remobilised in the subsurface layer

In the second situation recited above, the fraction of suspended particles may be computed as:

$F_{\mathrm{susp}}=\{\begin{array}{cc}\frac{{\tau }_{\mathrm{SUBSURFACE}}-{\tau }_{\mathrm{Cr},\mathrm{suspension}}}{{\tau }_{\mathrm{SUBSURFACE}}-{\tau }_{\mathrm{Cr},\mathrm{motion}}}& \mathrm{if}{\tau }_{\mathrm{Cr},\mathrm{suspension}}\ge {\tau }_{\mathrm{Cr},\mathrm{motion}}\\ \frac{{\tau }_{\mathrm{SUBSURFACE}}-{\tau }_{\mathrm{Cr},\mathrm{suspension}}}{{\tau }_{\mathrm{Cr},\mathrm{motion}}-{\tau }_{{\tau }_{\mathrm{Cr},\mathrm{suspension}}}}& \mathrm{if}{\tau }_{\mathrm{Cr},\mathrm{suspension}}<{\tau }_{\mathrm{Cr},\mathrm{motion}}\end{array}$

In the third situation recited above, part of the fraction of suspended particles is remobilized in the plume layer and may be computed as:

$F_{\mathrm{plume}}=\frac{H_{\mathrm{plume}}}{H_{\mathrm{plume}}+H_{\mathrm{sub}}}·F_{\mathrm{susp}}$

The remainder is remobilized in the subsurface layer:

$F_{\mathrm{subsurface}}=1-F_{\mathrm{plume}}$

With reference to FIGS. 8e to 8h, representing Shields diagrams expressing the shear stress values for remobilization of clastic and carbonates sediments according to their grain diameter, the critical motion shear stress depends on the grain diameter of the particle and can be determined using said Shields diagram. FIG. 8e relates to porous flat carbonates, FIG. 8f relates to porous spherical carbonates, FIG. 8g relates to non-porous spherical carbonates and FIG. 8h relates to clast. In these figures, the zone below the solid line is a zone where the particles are not remobilized, the zone between the solid line and the dotted line is a zone where the particles are remobilized by traction (i.e., in the bottom layer) and the zone above the dotted line corresponds to a remobilization by suspension. The solid line thus corresponds to a critical motion shear stress and the dotted line corresponds to a critical suspension shear stress. It should be underlined that for fine grains, for instance silt, fine sand or clay particles, there is no distinction between the critical suspension shear stress and the critical motion shear stress. As a consequence, these particles are either not remobilized or remobilized. If remobilized, these particles can be transported in the bottom layer if the shear stress value within the subsurface layer is lower than the critical suspension shear stress. Thus the equation for computing F_susp above is not applicable.

Last, for each layer, the volume of remobilized sediments during stormy events can be computed from the fraction of sediments remobilized in each layer (F_layer) multiplied by the volume of deposited particles in fair weather conditions V_{sed,depo,fairweather} (step 600):

$V_{\mathrm{sed},\mathrm{remobilized},\mathrm{layer}}=F_{\mathrm{layer}}·V_{\mathrm{sed},\mathrm{depo},\mathrm{fairweather}}$

In embodiments, and as schematically shown in FIG. 3b, the step 700 of remobilizing particles is followed by implementation of a step 600′ of modelling transport and/or deposition of the remobilized particle, which may also include a flocculation substep 601′. Said step 600′ is iterated until all particles are deposited or have exited the model.

Implementation of a step 600′ of modelling transport and/or deposition of particles after remobilization is performed the same way as a step 600 implemented prior to remobilization, except that the shear stresses computed for the modelled currents are computed in stormy conditions, i.e., taking {right arrow over (u_wind,storm)} as the wind speed vector for computing the velocity and shear stress values corresponding to the wind-induced currents.

During a step 800, the topography of the geological model of the area is updated to take into account the sediments deposited during the computational layer, i.e., an additional time layer is generated, which thickness is determined based on the volume of deposited sediments during the period of time corresponding to the computational layer. Step 800 may also comprise updating topography according to eustatism and subsidence, i.e., respectively the reference water level and the ground level may be changed by updating the position along z of each cell and the water depth of each cell.

At the end of the simulation, the modified model of the sedimentary area represents the evolved state of a sedimentary basin. This model may be used to extract data, such as for instance a log of sediments, to be compared with observed data acquired from a well.

The various embodiments described above can be combined to provide further embodiments. All of the patents, applications, and publications referred to in this specification and/or listed in the Application Data Sheet are incorporated herein by reference, in their entirety. Aspects of the embodiments can be modified, if necessary to employ concepts of the various patents, applications. and publications to provide yet further embodiments.

These and other changes can be made to the embodiments in light of the above-detailed description. In general, in the following claims, the terms used should not be construed to limit the claims to the specific embodiments disclosed in the specification and the claims, but should be construed to include all possible embodiments along with the full scope of equivalents to which such claims are entitled.

Claims

1. A computer-implemented method of modelling transport of a particle induced by water currents in an immersed area, wherein the immersed area comprises a plurality of cells associated with respective water depths, and the particle is located in a cell and represents a quantity of sediments of determined granulometry and sediment type,the method comprising: determining a direction and velocity of at least one water current occurring within the immersed area;determining, from the water current, a direction and intensity of a shear stress induced by the water current on the particle; anddetermining a transport of the particle, from the determined direction and intensity of the shear stress induced on the particle, the granulometry and sediment type of the particle.

2. The method according to claim 1, further comprising a preliminary step of defining three water layers corresponding to respective water depth ranges extending between a water surface and a water bottom of the immersed area, comprising: a bottom layer, located at the water bottom,a plume layer, located at the water surface, anda subsurface layer, extending between the bottom layer and the plume layer,wherein the method further comprises determining the direction and intensity of a shear stress in at least one water layer.

3. The method according to claim 2, wherein the direction and intensity of a shear stress in at least one water layer is determined based on a velocity profile according to depth of the at least one water current in the at least one water layer.

4. The method according to claim 2, further comprising determining a direction and intensity of a shear stress induced by gravity on particles located in the bottom layer.

5. The method according to claim 4, wherein the shear stress induced by gravity on particles is determined such that: the shear stress induced by gravity on the particle is null if a topographic slope of the water bottom is null;the value of the shear stress induced by gravity on the particle is equal to a deposition shear stress threshold when a current velocity in the subsurface layer is null and the topographic slope of the water bottom is superior or equal to an avalanche angle determined for the particle; andthe direction of the shear stress induced by gravity on the particle is parallel to a direction of a downward topographic slope of the water bottom.

6. The method according to claim 2, wherein the particle is located within one of the three water layers, and determining a transport of the particle comprises comparing the shear stress value in the water layer in which the particle is located to a shear stress threshold depending on the particle: when the particle is located in the plume layer or subsurface layer, and; when the shear stress value in the water layer in which the particle is located is higher than a suspension shear stress threshold, transporting the particle to an adjacent cell, said adjacent cell being determined based on the determined direction of the shear stress, orwhen the shear stress induced on the particle is lower than the suspension shear stress threshold, transporting the particle towards the subsurface or bottom layer, respectively, within the same cell; orwhen the particle is located in the bottom layer, and when the shear stress induced on the particle is higher than a traction shear stress threshold, transporting the particle to an adjacent cell, said adjacent cell being determined based on the determined direction of the shear stress within the bottom layer, orwhen the shear stress induced on the particle is lower than the traction shear stress threshold, depositing the particle.

7. The method according to claim 2, comprising a preliminary step of determining a thickness of an Ekman layer extending from the water surface of the immersed area, and defining the three water layers such that: the bottom layer extends between the water bottom and a fixed distance thereof,the plume layer extends between the water surface and a depth determined based on the Ekman layer's depth, andthe subsurface layer extends between the plume layer and the bottom layer.

8. The method according to claim 7, wherein the preliminary step further comprises defining parameters of wind occurring over the immersed area, and the method further comprises: determining the direction and velocity of a wind-induced current occurring in the plume layer based on the wind parameters; anddetermining the direction and velocity of a return current occurring in the subsurface layer and resulting from the wind-induced current occurring in the plume layer.

9. The method according to claim 8, wherein determining the direction and velocity of a wind-induced current occurring in the plume layer comprises determining the direction and velocity of an ocean surface current caused by an Ekman vortex.

10. The method according to claim 8, wherein determining the direction and velocity of a wind-induced current occurring in the plume layer comprises determining a direction and velocity of a wave-induced current.

11. The method according to claim 1, wherein determining the transport of a particle further comprises a preliminary modelling of aggregation of particles by flocculation, and the determining a transport of the particle is implemented on the aggregated particles.

12. The method according to claim 11, wherein modelling of aggregation of particles by flocculation comprises computing a proportion of aggregated particles based on a set of parameters including the determined shear stress and the particles size.

13. A computer-implemented method of modelling sedimentary deposition within an immersed area, comprising: a setup step, comprising:defining a geological gridded model of the immersed area,setting: a reference water level,at least one supply or production process of particles to be introduced within the model, andwind parameters, anda step of simulating an evolution of the geological gridded model over a predetermined period of time T, comprising: assigning a water depth to a plurality of cells;determining a direction and velocity of at least one water current occurring within the immersed area;determining, from the water current, a direction and intensity of a shear stress induced by the water current on a particle;introducing at least one particle in at least one cell of the geological gridded model;determining a transport of the particle, from the determined direction and intensity of the shear stress induced on the particle, the granulometry and sediment type of the particle; andupdating the geological gridded model of the area according to the transport of the particles.

14. The method according to claim 13, wherein the step of determining the transport of the particle is repeated until all introduced particles are deposited or have exited the geological gridded model.

15. The method according to claim 13, wherein each supply or production process is chosen among one of the following groups: clastic supply processes, comprising at least river mouth supply and mineral spring causing travertine deposition;carbonates production process; andeach supply or production process is associated with a depth at which a particle is introduced.

16. (canceled)

17. A non-transitory computer readable storage medium, having stored thereon a computer program comprising program instructions, the computer program being loadable into a processor and adapted to cause the processor to carry out, when the computer program is run by the processor, the method according to claim 1.

18. A computer, configured for implementing the method according to claim 1.

19. A computer, configured for implementing the method according to claim 13.