PAR field resolving to zero at depth despite surface PAR function definition

[TrialUser7]

Hi all and @jagoosw,

I am working on a 2D model (x,z) using Oceananigans, OceanBioME and SugarKelp. I am experiencing an issue where my SugarKelp particles are not photosynthesizing. I've diagnosed this because the kelp nitrogen reserves increase to the maximum limit (uptake is working), but carbon reserves decrease to the minimum (respiration is occurring, but photosynthesis is zero).

I have defined a surface PAR function using an interpolator from Excel data, but when I visualize the PAR field from my JLD2 output, the values are zero throughout the entire water column. I have ran the model with the default PAR and it works so I suspect there is a disconnect between my surface PAR definition and the light attenuation model.

I have tried troubleshooting but am wondering if there is a specific requirement for the PAR⁰ function signature or if the light field needs to be manually "forced" to update during the simulation to be used by the biogeochemistry?

Any help would be appreciated!

Below are some snippets of the relevant code from the model:

1. PAR Function and loading of data

# Load the PAR data
par_interpolator, par_times, par_values = load_par_data("PAR_data.xlsx")

@inline function PAR⁰(x, y, t)
    bounded_time = max(min(t, par_times[end]), par_times[1])
    return par_interpolator(bounded_time)
end

@inline PAR⁰(x, t) = PAR⁰(x, 0.0, t)
@inline PAR⁰(t) = PAR⁰(0.0, 0.0, `t)

2. Biogeochemistry

#Setting up Sugar Kelp parameters
kelp_model = SugarKelp()

particles = SugarKelpParticles(n; grid = grid,
                                advection = nothing, 
                                scalefactors = fill(scale_factor, n),
                                biogeochemistry = kelp_model
                                )
 set!(particles, A = initial_individual_frond_area, 
                N = initial_N, # gN/gSW
                C = initial_C, # gC/gSW
                x = x₀,
                z = z₀,
                )  

biogeochemistry = LOBSTER(; grid = grid,
                          light_attenuation = TwoBandPhotosyntheticallyActiveRadiation(; grid,
                          surface_PAR = PAR⁰),
                          surface_photosynthetically_active_radiation = PAR⁰,
                          particles = particles,
                          carbonate_system = CarbonateSystem(),
                          nutrients = NitrateAmmonia(),
                          detritus = VariableRedfieldDetritus(grid),
                          oxygen = Oxygen(),
                          scale_negatives = true,
                          sediment = nothing
                        )                                           

3. Output Files:

simulation.output_writers[:jld2] = JLD2Writer(model,
                                                    merge(model.velocities, model.tracers, 
                                                    (; PAR =model.biogeochemistry.light_attenuation.field)),
                                                    schedule = TimeInterval(0.1day), # Save every day
                                                    filename = "ocean_convection.jld2",
                                                    with_halos = false,
                                                    overwrite_existing = true)

simulation.output_writers[:particles] = JLD2Writer(model, (; particles),
                                                   filename = "kelp_particles.jld2",
                                                   schedule = TimeInterval(0.1day),
                                                   overwrite_existing = true) 

[jagoosw]

Hi @TrialUser7,

Thanks for reaching out, I'm trying to replicate your problem.

Is your grid Flat topology in y?

[jagoosw]

I tried it and I can't seem to replicate the issue (I set $PAR^0(args...) = 1$), what versions are you using? Perhaps this is an issue with your interpolator?