Smoothed particle hydrodynamics for modelling cold-water coral habitats in changing oceans

The importance of the growth, proliferation and longevity of reef-forming cold-water corals is paramount as they support various complex bio-diverse habitats and provide many essential ecosystem services. These cold-water coral reefs consist of layers of living coral tissue that grow on top of large masses of coral skeleton. Here, the Goldilocks Principle is used to simulate growth in optimal conditions and model how cold-water corals engineer their habitat to survive and prosper. A computational fluid dynamics model is created based on the Smoothed Particle Hydrodynamics method, a mesh-free Lagrangian numerical method. The SPH solver is written in the C++ programming language and parallelised with OpenMP to improve its efficiency and reduce the execution times. The solver is validated against analytical and numerical solutions and the growth model is then validated against in situ data of real cold-water coral colonies. The numerical results suggest that the longevity of cold-water corals depends on how well they can manage their energetic reserves when exposed to sub-optimal prey-catching conditions.


Introduction
Lophelia pertusa (see, Figure 1) is one of the most common species of framework forming cold-water corals (Roberts, 2006) that grows predominantly in the North Atlantic Ocean and has been found to form reefs world-wide.Typical L. pertusa reefs consist of live coral sitting on top of several layers of dead coral skeletons (Vad et al., 2017).Previous studies from (Pichon, 2011) and (Orejas et al., 2016) have demonstrated local flow hydrodynamics govern prey capture efficiency of coral tentacles; in low velocity environments food (including zooplankton) can evade capture, while in faster flow conditions coral tentacles are unable to capture food, as they are swept back by the flow (Purser et al., 2010;Orejas et al., 2016).Coldwater corals mainly satisfy their energetic demands by prey capture which for L. pertusa has been experimentally shown to be optimum when the local current velocity is between 2-6 cm/s (Tsounis et al., 2009;Purser et al., 2010).However, these corals typically exist in habitats with high current velocities that sometimes can be an order of magnitude higher than the experimentally found optimal velocity range.This leads to the question of how corals with such an optimal range can survive and thrive in the high flow conditions that they are found within.It has been assumed that cold-water corals build up lipid reserves during periods of high food availability (Dodds et al., 2009).They can then use these energetic reserves in periods that food availability is reduced, and Maier et al. (2019) demonstrated how L. pertusa can maintain their metabolic rate in periods of food deprivation.

J o u r n a l P r e -p r o o f
Journal Pre-proof Hennige et. al. (2021) explored the hypothesis that L. pertusa reefs are engineered according to the Goldilocks Principle (Figure 2).This assumes that coral polyps will survive and prosper if they surpass an energetic 'survival threshold' by capturing prey when conditions are optimal.This rule assumes that polyps can also survive in sub-optimal conditions if over time they capture enough prey to surpass the 'survival threshold' and cover their energetic demands.
Presently, an in-house developed Smoothed Particle Hydrodynamics (SPH) solver is used to evaluate coral growth based on the Goldilocks Principle.SPH is a mesh-free method, that uses particles to discretise the numerical domain.Traditionally in numerical simulations, a mesh of the domain must be created in order to create a discrete number of volumes, where mathematical governing equations that describe the physics of the flow can be solved to obtain the solution to a problem.This can create problems in simulations where the examined object is growing during the simulation.When this happens, the domain would have to be re-meshed, something that depending on the method can be complicated and time consuming, especially when mesh refinement algorithms need to be deployed.Conversely, this is not necessary in SPH, as all solids, fluids, and boundaries are represented by particles and the simulation can continue unaffected by changes in the boundary conditions due to evolution of an object (i.e., coral).SPH has its shortcomings too, as typically large resolutions are needed to capture a realistic representation of the flow field, which can increase execution times.For the coral growth model presented here, its benefits outweighed the disadvantages and is preferred to other mesh-based methods.
Figure 1: A picture of Lophelia pertusa framework illustrating live (white) and dead (grey) coral (Fox et al., 2016) J o u r n a l P r e -p r o o f

Journal Pre-proof
The growth and death rules that were firstly introduced in (Hennige et al., 2021) have been rewritten and optimised to include the effects of dynamic coral energetic reserves, while the parallelised SPH solver allowed for higher resolution simulations.The work presented here has two objectives:  Firstly, to validate the SPH solver against other numerical and analytical solutions.This is necessary to prove that it can solve the SPH governing equations accurately and provide information about the resolution that is needed to achieve the required accuracy. Secondly, to introduce the newly developed coral growth model, validate it, and examine different cases of coral growth.

Smoothed Particle Hydrodynamics
Smoothed Particle Hydrodynamics is a Lagrangian computational method that can be used to simulate the flow of viscous fluids.SPH guarantees conservation of mass without the need for extra computation (Liu and Liu, 2010).Its meshless nature allows natural tracking of fluidsolid interfaces.In the current work, this was fully exploited as the dynamically growing coral colonies imposed new boundary conditions after each growth step.SPH is an interpolation method so in order to update the properties of the particles the governing equations are expressed as summations of interpolants that use a kernel function, W, with smoothing length, h (Morris et al., 1997).Several conditions need to be met in order for a function to be considered appropriate as a kernel:  The function has to offer compact support, therefore: where k is a factor that defines and constrains the function's spread.This is necessary in order to reduce the computational cost of the kernel function. The function must meet the normalization condition: ∫  (   , ℎ)  ′ = 1  In order to avoid numerical instabilities, inaccuracies and unrealistic properties (for example negative density) the function has to be positive within its domain. The function has to offer symmetry, meaning that particles in equal distances to a reference particle should have the same contribution to its properties. Finally, the function has to ensure that the contribution of a particle to the properties of another particle reduces with increasing distances.
In the current work, the SPH solver uses a Wendland kernel function (Wendland, 1996) that can be more efficient than most cubic or quintic spline kernels (Macia et al., 2011) It was also shown that the dissipation mechanisms in Wendland kernels can be more accurate than those of re-normalised Gaussian kernels.The Wendland kernel function is defined as: where  = The SPH solver presented in (Hennige et al., 2021) was further optimised using OpenMP.The newly parallelised solver run up to 7 times faster and allowed for higher resolution simulations with reasonable execution times.Here, it solves the mass and momentum conservation equations; their discretised SPH forms are respectively: where ρi is the density of a particle i, mj is the mass of a neighbouring particle j, vij is the velocity difference between the two particles, ▽iWij is the gradient of the kernel function, pi is J o u r n a l P r e -p r o o f Journal Pre-proof the pressure of particle i, µ is the dynamic viscosity of the particles and Fi is an external force per unit mass.
Monaghan's method (Monaghan, 1994) of approximating the rate of change of density is being used in the current work for the computation of the particles' density.According to this method the particles are initially set to a reference value and their density evolves by solving the continuity equation (3).After the density computations, a density correction algorithm is applied (Ozbulut et al., 2014).In weakly-compressible SPH the pressure of particle is calculated using an artificial equation of state and it is directly connected to the particle's density.Therefore, a density correction algorithm helps to avoid large density variations in the domain that can lead to numerical instabilities and inaccuracies.The density is being smoothed using: ( 5) The current work includes relatively small velocities, and the particles fill all available space, therefore a realistic form of viscosity was adopted as suggested by Morris (1997) and seen in the Navier-Stokes momentum equation (4).
In SPH, pressure is a function of density, and the movement of the particles is driven by density fluctuations and consequently an artificial equation of state has to be used.The equation of state for water (Ree, 1976) could be used as well, but that would require incredibly small timesteps (Morris et al.,1997) making the simulations inefficient.Here, the Tait equation was used: The value of the polytrophic constant γ must be chosen carefully in order to ensure the accuracy of the solution; for water, γ = 7 The initial pressure,  0 , depends on the reference speed of sound for the fluid according to: In this work, suggestions by Monaghan (1994) and Violeau (2000) have been used to define the speed of sound, c, which should be at least 10 times greater than the maximum velocity in the domain.This can reduce density fluctuations in the domain to within 1% of the reference density of a particle (Monaghan, 1994).
A particle shifting algorithm is used to avoid stability issues caused by anisotropic particle spacing.This algorithm moves particles to areas with lower particle concentration in order to avoid the creation of voids and maintain a uniform distribution throughout the domain.Here, the algorithm proposed by Skillen and Lind (Skillen et al., 2013) is used.In this algorithm, the shifting distance, δr, is given by: J o u r n a l P r e -p r o o f Journal Pre-proof where C is a concentration coefficient and D is a diffusion coefficient and that can be calculated by: where dt is the time-step of the simulations, |v|i is the velocity magnitude of a fluid particle and h is the smoothing length.Finally, the gradient of the concentration coefficient in equation ( 8) can be calculated using: The algorithm can struggle in simulations with free surfaces, where a correction is needed (Skillen, 2013) but since the current work involves only internal flows this is unnecessary.
A Verlet scheme (Verlet, 1967) coupled with a Euler step (Jameson et al., 1981) every 50 iterations has been used to perform time integration.The Euler step is necessary to ensure that the equations remain coupled for odd and even time-steps and time-integration divergence is avoided.In order to ensure the stability of the simulations, the time-step is calculated using the Courant-Fredrichs-Lewy (CFL) condition (Liu and Liu, 2010) and two additional restrictions to account for viscous dissipation and body forces (Monaghan, 1994).
The seabed and coral solid surfaces are simulated using dynamic boundary particles (Crespo, 2007).The positions and velocities of these particles remain fixed over time.The motion in the numerical domain is driven by the moving upper boundary that consists of three layers of dynamic boundary particles with their velocity fixed at 0.5 m/s to simulate the typically fastflow environment that cold-water corals grow in.In total, 80,000 particles are used with initial particle spacing equal to ∆= 0.025m.

Basic Growth Principles
This novel long time-scale growth model was created in order to investigate how cold-water corals would grow according to the Goldilocks Principle.During a growth cycle, the average local steady-state flow velocities of fluid particles that are in close proximity (∆ < 1.5, where r is the initial particle distance) to any coral particle are analysed.If the velocity of any fluid particle adjoining a coral particle lay inside the Goldilocks zone, then that fluid particle is converted into a coral particle.No additional particles are inserted or deleted from the numerical domain.Essentially, the model is looking for zones of optimal velocity around a coral colony to identify the direction of growth of the colony.No additional rule is applied to control branching; it occurs spontaneously where the flow is optimal.
The introduction of a death rule in the model was a significant step to simulating coral growth.
In previous work (Hennige, 2021), the death rule was fixed at initialisation, and it could not be altered or affected by any other aspect of the live simulations.The current proposed model alters the way that the growth and death rules affect the coral particles.

J o u r n a l P r e -p r o o f
Journal Pre-proof Each coral particle's energetic reserves receive a finite value of energy units when initialised.
One unit of energy is assumed to be equal to the energy that a coral particle would need to survive between two growth steps in sub-optimal conditions.The energetic reserves decrease or increase dynamically according to the flow conditions around the coral particles.For example, in sub-optimal conditions a coral particle is not able to 'catch enough prey' and get the energy it needs; therefore, it has to use a portion of its energetic reserves in order to survive.
If the conditions around the coral particle do not improve over time and the value of its energetic reserves drops below zero, then the coral particle was considered to be dead.This can be shown in the following equations: where ER is the energetic reserves, ts the current time step and ts-1 the previous time step.
Conversely, if the flow conditions are right then the coral particles are able to get enough energy to survive and grow and could potentially be able to increase their energetic reserves according to: The quantity θ ranges between 0 and 1 units of energy.A value θ equal to 0 units of energy means that all energy that is generated by capturing prey between two growth steps is used by the coral particle to satisfy its energetic demands and no additional energy can be stored.
Similarly, a value of 1 means that all generated energy can been stored to the energetic reserves of the coral particle.

SPH Coupling with coral growth model
A typical SPH coral growth simulation would involve the following:  Initially the geometry, boundary conditions, and input conditions are provided. A particle neighbour list is obtained, for every fluid, boundary, or coral particle in the domain. The solver can then solve the Navier-Stokes equations and update the properties of all particles, as explained in Section 2.1 above. Next is an important step for coupling the SPH solver with the coral growth model.The solver will determine whether the flow in the numerical domain is in steady-state conditions or not.If this is true, it will proceed to instigate the coral growth and death functions, as explained in Section 2.2 above.The solver will determine flow conditions around live coral particles, and where appropriate it will simulate coral growth towards directions of optimal flow velocity.Additionally, the energetic reserves of every coral particle will be re-calculated and if a particle in sub-optimal flow conditions has no energetic reserves left will be considered dead.

J o u r n a l P r e -p r o o f
Journal Pre-proof  The modelled coral will grow, providing thus new boundary conditions for the numerical domain. The solver will then proceed to the next time-step, obtaining a new particle neighbour list, solving the Navier-Stokes equations, and when the flow is again in steady-state conditions, the growth functions will once again be instigated.
This procedure can also be seen in the following algorithm in the form of pseudo-code.

J o u r n a l P r e -p r o o f
Journal Pre-proof

Results and Discussion
Before presenting the results of the coral growth model, a few validation cases are provided to ensure that the presented methodology can provide accurate solutions when compared against other known numerical or analytical methods.In all of them, particle convergence tests were performed to identify the needed resolution to achieve the required accuracy.

Poiseuille Flow
The first validation case considered a Poiseuille flow problem, where the flow in the domain was driven by a pressure gradient force.The fluid (water) in the domain was placed between two stationary plates with infinite length.The testing setup proposed by Morris (1997) was used for this case, the properties and initial conditions are shown in Table 1.The two stationary plates were initialised with a distance of 0.001m between them and consisted of three layers of dynamic boundary particles.For this problem a Wendland kernel was used, and the speed of sound was chosen to be 100 times larger than the maximum velocity in the domain.Periodic boundary conditions have been applied on the left and right boundaries in order to simulate an infinite domain.The initial separation between the SPH particles was dx=dy and it depended on how many particles were spanning the channel between the two stationary particles.It can be calculated according to: where L is the distance between the two plates and N is the number of particles in the y direction.
The analytical solution for the Poiseuille flow was obtained by using the equation (Morris, 1997): where   is the velocity of the water in the x-axis, ν is the kinematic viscosity of water, ρ is the density, t is the elapsed time and n is the number of terms to include in the summation.

J o u r n a l P r e -p r o o f
Journal Pre-proof Figure 3 compares the numerical and analytical solutions.Good agreement was found between them, with the maximum error in the numerical solution being close to 1.1% when 100 particles were spanning the channel.The error in the simulations with various number of particles spanning the channel between the two stationary plates can be seen in Table 2.As it can be seen 100 particles in the y-direction were enough to achieve particle convergence as the error in the simulations did not decrease significantly after that, while the computational cost was increasing with more particles in the domain.

Couette Flow
The next validation case was a two-dimensional Couette flow, which is a flow between two infinitely long plates where the bottom plate is stationary, while the top plate has constant velocity.It was an interesting case for a system where viscous dissipation was important, a similar system to what the coral model simulations would use.The test case was developed using the setup proposed by Morris (1997).The initial conditions are shown in Table 3 below.The two plates where again placed at distance equal to L = 0.001m apart, but in this case only the bottom plate was stationary.Throughout the simulation the top plate had constant velocity equal to 1.25x10 -5 m/s.Both plates were simulated using three layers of dynamic boundary particles.The analytical solution for the simulated two-dimensional Couette flow was obtained by using the equation: where  0 is the velocity of the moving top plate.
Figure 4 shows a comparison with the analytical solution.The maximum error in the numerical solution was less than 0.7% in simulations with 100 particles spanning the channel between the plates, showing good agreement with the theoretical results.A particle convergence test was conducted again and showed that 100 particles in the y-direction were enough to consider that the simulations had converged (Table 4).

Lid-driven Cavity Flow
The next validation case showed the capability of SPH to model flows in higher Reynolds numbers.This case simulated the flow inside a square cavity that had its lid (top boundary) moving with constant velocity,  0 .The fluid inside the cavity was initially at rest and it started moving due to viscous forces caused by the movement of the lid.The square cavity consisted of four walls of equal length, L, and each wall had three layers of dynamic boundary particles.
A Wendland kernel was used and the initial distance between the particles depended on the resolution of the simulations according to:

J o u r n a l P r e -p r o o f
Journal Pre-proof where N is the number of fluid particles per direction.For this case, the number of particles per direction ranged from 50 to 220 particles.
The speed of sound was chosen to be 100 times larger than the maximum velocity in the system (the lid's velocity).A schematic of the case can be seen in Figure 5 and all initial properties and parameters of the simulations can be seen at Table 5 below.Two different cases were run for Reynolds number equal to Re = 1000 and Re = 10000.The Reynolds number in the simulations was adjusted by modifying the viscosity of the fluid, while the density of the fluid, the characteristic length (L) and the maximum velocity were kept constant.Ghia (1982) who used a finite volume solver with a 257x257 mesh and results by Adami (2013) who used a weaklycompressible SPH solver with transport velocity formulation.There is no analytical solution available for this case.

Property Units Value
The results for Re = 1000 can be seen in Figures 6-8, while the velocity field obtained by Adami ( 2013) is shown in Figure 9. Similarly, Figures 10-12 show the obtained results for Re = 10000 and Figure 13 shows the corresponding velocity field by Adami.The velocity fields shown in Figures 8 and 12 are at time-steps that the cases had reached steady-state conditions.Adami did not provide a legend for their velocity fields, but it can be assumed that it is similar as the one in Figures 8 and 12.It was found that in low resolutions (N<150) the current SPH solver struggled and required higher resolutions in order to provide meaningful comparisons.Therefore, only results from higher resolution simulations are shown.The results of the SPH solver showed good agreement with results obtained by the other solvers.For Re = 1000 the necessary accuracy was achieved by using 220x220 particles in the domain.The same resolution was necessary for Re = 10000, but in this case additional boundary particle treatment had to be performed in order to ensure that the particles will not escape the numerical domain.This also increased the accuracy of the solution, but at the cost of additional execution time.
An additional repulsive force between the fluid and boundary particles was added, as suggested by Monaghan (1988).For this work, this was achieved by adding the following Lennard-Jones force term in the Navier-Stokes momentum equation ( 4): 18) where  1 = 12 and  2 = 6 are constants, dx is the initial particle separation and D is equal to 120 times the product of the initial particle separation and the acceleration due to gravity.
Equation ( 18) prevents fluid particles from penetrating the solid walls.

Coral Growth Model
The growth model investigated how energetic reserves can affect coral growth and mortality.
The simulation parameters (Table 6) are based on a mono-directional flow from left to right and a simplified growth principle existed; the coral colony would only grow in optimal conditions and towards regions with average flow velocities between 2-6 cm/s.It investigated and showcased how the Goldilocks Principle can be applied to cold-water coral growth and how coral energetic reserves can affect their growth and longevity.These cases were run with the value of θ in Equation ( 13) being kept constant at 0.5.∞ 1.1 3.1 Ratio of live to total coral particles, (at step 120) 100 0 10.7 Initially, (Figure 14, (A)) a control case was simulated; the model only simulated growth in optimal flow conditions with no additional `death' rule applied.In sub-optimal regions the coral would not grow but also not die, therefore simulating infinite energetic reserves.The modelled coral in this simulation would grow indefinitely and cover the simulated domain.This was a direct contradiction to what can be observed in nature where a significant portion of L. pertusa reefs consists of calcified dead coral skeleton (Roberts et al., 2009).
Previous studies (Larsson et al., 2013;Baussant et al., 2017) have shown that L. pertusa reefs can survive in sub-optimal conditions for a period of months by using their energetic reserves to cover their energetic demands.The next model (Figure 14, (B)) introduced a `death' rule that was based on each coral particle's available energetic reserves (Table 6).The amount of the initial energetic reserves for a coral particle indicated how many consecutive growth steps this particle could survive in sub-optimal conditions.If during a growth step the coral particle was in sub-optimal regions and had no available energetic reserves left, then the death rule was applied, and the particle was considered 'dead'.It could no longer branch out and grow but it would be part of the coral's skeleton for the remaining simulated time.When the coral particles were initialized to have near-zero or very low energetic reserves (Figure 14, (B)) then eventually the entire coral framework died when it was exposed to consecutive non-optimal flow conditions.In the simulations with higher initial energetic reserves (Figure 14, (C)) the

J o u r n a l P r e -p r o o f
Journal Pre-proof resulting coral framework consisted of dead coral skeleton on the inside with branches of live coral on the outer edges, similar to what can be observed in real L. pertusa reefs (Figure 1).

Replenishing energetic reserves
The main equation that controlled the growth and death rules in the models (Equation 13) also offered the capability of running simulations where the energetic reserves were not static and predetermined, but dynamically altering for each individual coral particle.The value of θ in Equation ( 13) controlled the proportion of the energy created during optimal time-steps, that was stored to the energetic reserves of the coral particles.For example, a value of θ = 0.5 would mean that a coral particle in optimal conditions would use half of the energy it obtained by catching prey to meet its energetic demands and store the other half.In the models presented in this subsection, coral energetic reserves were tracked individually for each live coral particle in the domain.
The growth model was run again with increasing abilities of replenishing energetic reserves on optimal steps, θ, in order to investigate its effects on the longevity of the coral colonies.The results are presented on Table 7 below.It demonstrates how the ratio of live coral particles to the total coral particles in the domain was affected by the increasing ability of the particles to replenish their energetic reserves.The relative average energetic reserves of the live coral particles are shown as well.This is calculated as the average energetic reserves of all live coral particles at that specific growth-step divided by the initial energetic reserves (ER in equations ( 11)-( 13)).This ratio was used to enable direct comparison between simulations that were initialized with various values of initial energetic reserves (ER).12.5 ± 0.13 1.78 ±0.09 0.9 14.1 ± 0.14 1.96 ±0.09As expected, when coral particles had higher abilities to replenish their energetic reserves (higher values of θ) the resultant coral colonies had higher average energetic reserves.It is also notable that being able to stay alive for longer resulted to colonies that had higher number of live coral particles compared to the total amount of coral particles in the domain.Table 7 shows that this ratio was dropping as the simulations progressed and was higher for simulations that allowed the colonies to replenish their energetic reserves faster (simulations with higher θ values).
Cold-water corals are characterised by various processes that require high energetic inputs; calcification, tissue and mucus production, reproduction and maintenance (Hennige et al., 2014).In more acidic conditions, the energetic demands associated with calcification rates could be higher, with more energetic reserves used to maintain stable calcification rates (Hennige et al., 2014(Hennige et al., , 2015)).The SPH model presented in Table 7 examined how the rate that L. pertusa can replenish energy during growth steps with optimal flow conditions can affect J o u r n a l P r e -p r o o f Journal Pre-proof coral longevity.The results suggest that when the coral particles were allowed to replenish more energy in optimal time-steps (higher θ values), colonies had a higher ratio of live coral particles to total coral particles in the domain.As expected, this also meant that in higher θ value simulations the average energetic reserves at later stages of simulations were higher and these colonies could therefore survive longer in sub-optimal conditions.The long-term prosperity and longevity of the coral colonies therefore would depend on their ability to store a portion of the energy they create by capturing prey.This could make a significant difference in periods that they are exposed to continuous sub-optimal flow conditions or in situations where their energetic demands increase due to changes in environmental variables such as in more acidic waters (Secretariat of the Convention on Biological Diversity, 2014; Hennige et al., 2015).12.5 ±0.13 0.9 16.9 ±0.16 14.1 ±0.14A previous study (Vad et al., 2017) examined various L. pertusa colonies from two different sites and showed that the ratio of living coral to the whole colony size was between 0.10 and 0.27.It was also shown that this ratio is negatively correlated to the whole colony size.Table 8 shows the ratio of live coral particles to the total coral particles in the domain in simulations with various energetic reserve configurations.In the modelled coral colonies, the ratio varied between 0.098 and 0.17 between the 50th and the 100th growth step showing the same negative correlation with the colony size as well -the ratio drops in value as the simulations progress while the total coral particle number can only increase.
At later stages of the simulations the ratio of live coral particles to the whole coral particles is lower than what observed by (Vad et al., 2017).The domain size was chosen initially to be large enough that it would not affect the growth of the coral colonies, but also not too large that it would make the simulations very computationally expensive.As the coral colonies grow and occupy larger parts of the finite numerical domain, it is possible that at later stages the dimensions of the domain start to affect growth.The ratio of live coral to total coral particles can be used then as a method to end the simulations, when it starts to drop too far below the expected values.

Coral growth and gravity
A limitation of the previous SPH model (presented in Figure 14) was that it disregarded the effects of gravity in coral growth.The coral would grow in the nearby optimal velocity regions and new layers of live particles would be created on top of previous layers.In real colonies it would be impossible for a single point to support all this newly created mass of coral structure above it and the colonies would start to break-down according to the mechanism suggested by (Wilson, 1979).In order to mimic this mechanism and visualize more realistic coral colonies J o u r n a l P r e -p r o o f an intermediate 'gravity' step has been included in the following model.Here, the coral colony would initially grow similarly to the previous model until it reached the 60th growth step.At that moment a break-down mechanism was initiated to simulate the effects of gravity to the coral colony.After the simulation of this intermediate gravity-step reached steady state, the additional gravitational acceleration was again set to zero, the particles of the top boundary were re-initialised with the input velocity (0.5 m/s) and the growth model started once again to simulate coral growth based on the newly imposed boundary conditions.The results can be seen in Figure 15 below.
Figure 17 shows the velocity vector shortly after the gravity step of the growth model (70th growth step) while Figure 18 visualizes the velocity profile near the end stages of the simulation (110th growth step).The velocity magnitude of the water in the domain is zero at the bottom boundary where the noslip condition is enforced, and it increases with the height of the domain until it reaches its maximum value (0.5 m/s) at the top boundary (omitted in the figures).A region of recirculating   The survivability and longevity of cold-water coral colonies depend on how they manage their energetic reserves.In growth-steps where they are exposed to sub-optimal flow, they need to have enough energy stored to allow for the smaller inflow of resources to prevent mortality.In growth-steps within optimal flow regions they need to store enough energy to ensure that their reserves are not depleted, and they can survive potential future sequential growth-steps in suboptimal conditions.This highlights the importance of coral energetic reserves; they are shown to be one of the major factors that can affect coral growth and prosperity.Management of energetic reserves is paramount in periods where their energy intake is decreased, or their energetic demands are increased due to changes in environmental variables.
The outputs of the models are in accordance with in situ studies that compare the size of the living coral in a colony to the size of the whole colony.The modelled corals show similar growth patterns as real cold-water corals and the ratio of living coral to the total colony size is negatively correlated with the size of cold-water coral colonies.Furthermore, qualitative comparisons against real cold-water coral colonies illustrate similar dense, complex branching geometries with high rugosity at the outer edges.They consist of a layer of living coral particles surrounding dead coral skeleton.
The Growth model in this work considered only the effects of hydrodynamics in cold-water coral growth, assuming that the available nutrients are infinite.A more realistic approach would be needed in order to capture the effects that nutrient availability can have in coral growth, where upstream nutrient uptake can affect downstream availability.Decreasing food availability could lead to less symmetrical coral growth forms with upstream positions having an inherit resource advantage.This would open the way to model competition among multiple coral colonies for the finite resources.The development of this SPH model can also lead to modelling various future scenarios, including the effects of ocean acidification on coral framework and potential coral restoration practices.Additionally, the methodology can be extended to model tropical coral growth by introducing sunlight as an input growth variable.

Competing Interest
The authors declare that they have no competing interests.

Figure 2 :
Figure 2: The Goldilocks Principle for Lophelia pertusa, showing cumulative prey capture over time, compared to local current velocity conditions.The bisecting layer indicates a 'survival threshold' that is based upon prey capture.Not surpassing this threshold would lead to polyp mortality, leaving exposed 'dead' framework.Individual polyps can surpass the threshold by either capturing prey in optimal conditions or sub-optimal conditions given adequate time.J o u r n a l P r e -p r o o f

Figure 3 :
Figure 3: Comparison between the SPH numerical solution (spheres) and the analytical solution obtained with Equation 15 (solid lines) at different time-steps J o u r n a l P r e -p r o o f

Figure 5 :Figure 4 :
Figure 5: A schematic of the lid-driven cavity flow, showing the moving lid (red) and the solid stationary walls (black).The length (L) of each side of the square is equal to 1m

Figure 6 :
Figure 6: Velocity profile Vy(x) at the centre-line y = 0.5m, Re = 1000 J o u r n a l P r e -p r o o f

Figure 9 :
Figure 9: Velocity field (only fluid particles), Re =1000 Figure 8: Velocity field by Adami (2013) J o u r n a l P r e -p r o o f

Figure 14 :
Figure 14: Coral Growth using the Goldilocks Principle in simulations with A) infinite initial energetic reserves, B) low initial energetic reserves and C) high initial energetic reserves.Growth is shown between the 20th, 60th and 120th growth steps.

Figure 15 :
Figure 15: Including the effects of gravity in coral growth.Initially a coral colony grew from a single point for 60 growth steps (A).At this point a break-down mechanism was initiated and the resultant colony is shown to include gravitational effects (B).Finally, growth in the domain has been re-initiated and coral growth at 120 growth steps is shown (C).

Figure 17 :
Figure 17: Velocity profile with streamlines around the coral colony at the 110th growth step.Red particles show live coral particles while yellow particles denote dead coral framework

Figure 16 :
Figure 16: Velocity profile with streamlines around the coral colony at the 70th growth step.Red particles show live coral particles while yellow particles denote dead coral framework.J o u r n a l P r e -p r o o f This work was supported by a Natural Environment Research Council (NERC) Doctoral Training Partnership (grant no.NE/L002558/1) to K.G, and Operation Wallacea.This work was supported by the Independent Research Fellowship from the Natural Environment J o u r n a l P r e -p r o o f

Table 1 :
Initial properties of the SPH particles in the Poiseuille flow validation case

Table 2 :
Particle convergence test for the Poiseuille flow, showing the number of particles spanning the channel between the two stationary plates and the corresponding error between the numerical and analytical solutions Number

Table 3 :
Initial properties of the SPH particles in the Couette flow validation case

Table 4 :
Particle convergence test for the Couette flow, showing the number of particles spanning the channel between the two stationary plates and the corresponding error between the numerical and analytical solutions

Table 6 :
Coral growth simulation parameters

Table 7 :
Dynamic energetic reserves in 2D simulations.The presented properties of the coral colonies are taken from the 100th growth-steps of the simulations

Table 8 :
Ratio of live to total coral particles in the domain at the 50th and 100th growth steps based on the simulated ability of the colonies to replenish their energetic reserves