Resolving issues with expanded number of coral species

[ConnectedSystems]

Background

@KRNA01 and I have been plugging away at expanding the number of corals represented in ADRIA.

In the previous version, there were four coral "types"

In this updated version, there are now 36 coral "types", consisting of 4 taxonomic groups + 2 enhanced, which are in turn split into six size classes:

• Tabular Acropora Enhanced
• Tabular Acropora Unenhanced
• Corymbose Acropora Enhanced
• Corymbose Acropora Unenhanced
• Small Massives
• Large Massives

I have managed to vectorize the growth function such that loops over sites are no longer necessary. Only the loop over time steps remain.

The Issue

Running the coral growth function through ode45() eventually returns Inf (usually around time step 10 - 15), which leads to NaNs in the subsequent time step, rendering all subsequent results unusable.

If I am reading the code correctly, the function takes coral populations at time t

• factors in recruitment for the smallest size classes, and
• competition between Tabular Acropora and Small Massives

... to arrive at a population size for time t+1.

Below are scenario results when the function is used directly. To me, this is in line with expected model behaviour (at least more so than infinite growth), although admittedly I have a very naive view of what "expected" is. Due to these results below, I'm now not so convinced that use of ode45() is necessary, at least in this specific context.

I think if it were only the growth function then we could use ode45() in place of the time step loop. We do, however, have the decision making procedure to include so I'm not sure its possible.

@KRNA01 , I think the other earlier concern about interactions between species not being represented when looping over species and using the function directly as with my earlier attempt is now resolved as:

• we are no longer looping over species
• inputs to the functions are for all species at time t so the interactions between species (e.g., competition) are being represented

I think the confusion was due to poor explanation on my part on how those function handles work, sorry.

Useful links in case I am entirely wrong and use of ode45 is necessary:

• https://au.mathworks.com/matlabcentral/answers/519759-why-ode45-returns-nan-values

[KRNA01]

@ConnectedSystems, @Rosejoycrocker and @veroniquelago we have now a set of alternative ODEs for testing. The one I have most confidence in is growthODE4().

Below are the results for the 'coral_example' in the sandbox using the 8 default intervention combinations.

Without checking those intervention combinations, the behaviour appears in line with what I'd expect. Will dig deeper into this in the coming days.

Please note that the recruitment function has yet to be updated. Here, we need to calculate larval productivity for each of the six coral 'species' (four taxa) across size classes. This will likely become an extension of 'larval_prod()', if you agree.

Also, once this works, we can replace the currently ecosystem-service (ES) translations with the ones we have developed for ReefMod based on the expert data. This means our projected scope-for-ES provision will then speak more directly to tourism, fishing and non-use (existence benefits).

[ConnectedSystems]

Hi @KRNA01

This looks awesome, thanks. I've noted that you've adjusted the default RCP scenario value to "45", but in the older ADRIAparms file. I've pushed changes to match this in the more recent simConstants file. Please note this file is also a precursor to moving everything to the *Details function format I outlined elsewhere previously - still thinking on how to support parameter sampling.

Also, I've disabled the debug check/plotting. The entire example now runs in ~2.3 seconds (0.1 seconds per scenario) on a single core.

[KRNA01]

@ConnectedSystems, have now converted Seed1 and Seed2 from number of corals seeded to their relative coral cover (proportions). Not sure I have done this in the right place.

The conversion is simply number of corals seeded times their planar surface area (pi *1.75 cm radius squared) divided by 10^4 cm2 per m2 divided by the 100 m2 reef arena.

[KRNA01]

@ConnectedSystems, correction:

The conversion is simply number of corals seeded times their planar surface area (pi *1.5 cm radius squared) divided by 10^4 cm2 per m2 divided by the 100 m2 reef arena.

[KRNA01]

@ConnectedSystems and @Rosejoycrocker , now one step closer to what I think is a working version of the new coral function. At some point soon, I'd be grateful for a rigorous review of the steps, especially of 'coralDetails()' and 'growthODE4()' .

Below are the results of a 'single_scenario()' that @ConnectedSystems set up (very useful- thank you!). Behaviour is in line with what you'd expect, except that baseline (counterfactual) trajectories are more pessimistic (using published ReefMod parameter values for growth and mortality) than I would have expected. The run below is for ten years of coral seeding (1000 per 100 m2) for enhanced Tabular and Corymbose Acropora and ten years of shading under RCP 4.5/ We still need to update 'larval_pro()' for the expanded coral function, and this will add to more positive trends to the unenhanced corals. Trajectories are for the 26 sites. You'll see that less than 10 sites were selected for enhancement.

[ConnectedSystems]

FYI @KRNA01

After some minor changes to growthODE4() for performance

Single core:

With parfor:

Incidentally, it would be good to come up with more indicative names for the growthODE group of functions or at least document what the differences are in the docstring. Lets discuss in the new year :)

[KRNA01]

Hi @ConnectedSystems and @Rosejoycrocker,

Do you agree that we have solved this issue? If so, I suggest we close.

[ConnectedSystems]

Agree, I believe all work has now been completed on this and merged into 40corals.

Closing 🥳