Adds a mixed layer depth diagnostic

[jagoosw]

Hi @tomchor,

This PR adds a simple mixed layer depth diagnostic using a simple density anomaly criteria.

I don't think there is a convenient way to generally write a mixed layer depth computation, so I have called it DensityCriteriaMixedLayerDepth.

[tomchor]

Nice! I'm okay with the name, although I think it should be density criterion, and not criteria (single it's one, singular, criterion, no?).

Another alternative is to create an object called MixedLayerDepth that takes a criterion (or algorithm to be more precise) as the computation. That would make it easier in the future to have other way of calculating the MLD. For now the only algorithm would be the density threshold criterion.

Also, I see you're using a KernelFunctionOperation, but you want to return a 2D object, no? Is this the best way to do things? (I'm genuinely asking; I'm not sure.)

[tomchor]

(Also github asked for explicit approval to run the tests and docs. I have no idea why, but we should figure out how to disable that!)

[jagoosw]

Nice! I'm okay with the name, although I think it should be density criterion, and not criteria (single it's one, singular, criterion, no?).

Ah yes you are right, I have changed now.

Another alternative is to create an object called MixedLayerDepth that takes a criterion (or algorithm to be more precise) as the computation. That would make it easier in the future to have other way of calculating the MLD. For now the only algorithm would be the density threshold criterion.

I've done this now, I'm not sure which is better but you can have a look and I can roll back the commit if you prefer the first way.

Also, I see you're using a KernelFunctionOperation, but you want to return a 2D object, no? Is this the best way to do things? (I'm genuinely asking; I'm not sure.)

This is a good question, it seems to work for 2D objects by specifying a nothing location like KernelFunctionOperation{Center, Center, Nothing} and you can index with just two indices like mld[i, j]

[jagoosw]

(Also github asked for explicit approval to run the tests and docs. I have no idea why, but we should figure out how to disable that!)

I think this is the default behaviour for first-time contributors to a repo incase the CI has some kind of cost associated with it

[tomchor]

@jagoosw Sorry, I realize I didn't check this PR after your latest pushes. Conflicts aside (I can resolve them for you), is this PR close/ready to merge?

[jagoosw]

No problem, I forgot all about this too! I think its all ready if you're happy with it apart from the conflicts

[tomchor]

No problem, I forgot all about this too! I think its all ready if you're happy with it apart from the conflicts

Nice! Good to hear. I'll work on the conflicts today and take a closer look. If all is good, I'll approve and then you're free to merge. I invited you to become a collaborator here if you want, that way you can merge your own PRs.

[tomchor]

Ah, I just realized you're using your own fork to implement this so I can't modify it!

I think all you need to do is to move the test from runtests.jl to test_tracer_diagnostics.jl for now and then we can go from there.

[tomchor]

@jagoosw if you'd like me to fix things you can make me a collaborator on your fork I guess? Not sure if there's a better way since this PR is already created based on it.

[jagoosw]

I'll add you to my fork and try and clean these suggestions up

[glwagner]

You should not use a density criteria as it is not general. A buoyancy criteria is better.

[tomchor]

The way I understand the code, different criteria can be plugged into MixedLayerDepth, no? Having a density criterion doesn't preclude having a buoyancy criterion.

That said, can you elaborate why density criteria are less general than buoyancy criteria?

[glwagner]

That said, can you elaborate why density criteria are less general than buoyancy criteria?

Gravitational forces enters a Boussinesq model via buoyancy. One can simulate buoyancy directly, leading to the formation of "mixed layers". But density is not defined in such a model, so density-based criteria cannot be used.

More generally, often mixed layer depth criterion actually use some potential density rather than in-situ density, because they do not want to incorporate the dynamically irrelevant part of density due to compressibility.

Using buoyancy is simpler and avoids those issues; by design, the buoyancy field can be computed from any model that has gravitational forces.

It's nice to have different criteria but I would discourage developing any criteria other than ones based on the buoyancy profile.

[tomchor]

@jagoosw I moved the test to its proper place, but it's failing with some version of

Tracer diagnostics tests: Test Failed at /home/tomas/repos/Oceanostics.jl/test/test_tracer_diagnostics.jl:194
  Expression: mld[1, 1] ≈ -zₘₓₗ
   Evaluated: 1.3877787807814457e-17 ≈ -0.5

Can you please take a look?

[jagoosw]

I've added a buoyancy criterion now, and it makes sense that we should prefer that since it will work in all models, whereas the density one only worked with BoussinesqEquationOfState.

The test also passes locally for me now, but I have no idea what made it work again. I added a test for the buoyancy criterion and made it the default also. I think it would be nice to keep the density criterion for comparison to observational work but am not attached to the idea.

[glwagner]

The density anomaly can be expressed in terms of a buoyancy anomaly, where you require additional information about 1) the gravitational acceleration used in the model and 2) the reference density. For example:

struct DensityAnomalyCriterion
    buoyancy_anomaly
    reference_density
    gravitational_acceleration
end

You can then build a constructor that uses some default values for the parameters.

This will help users understand how the parameters they are using affect their results leading to higher-quality science.

[jagoosw]

I like this more

[glwagner]

Suggested change

	pertubation :: FT = - 1/8 * 1020 / g_Earth
	anomaly :: FT = - 1/8 * 1020 / g_Earth

[tomchor]

The test also passes locally for me now, but I have no idea what made it work again.

It seems like the errors are still there! I wonder why you can't reproduce these errors locally.

[jagoosw]

hmm when I run the whole file it does fail too now, it looks like its the stretched grid cases, I will investigate now

[jagoosw]

Ah its because the stretched grid doesn't have z=0 at the top

[tomchor]

I like the science part here. However, I'm a little confused by the implementation. MixedLayerDepth is defined as a a struct which takes only a criterion. But then if one calls something like MixedLayerDepth(grid, buoyancy, C; criterion), then the following method is called:

Oceanostics.jl/src/FlowDiagnostics.jl

Lines 467 to 473 in f9ac0e1

	function MixedLayerDepth(grid, args...; criterion = BuoyancyAnomalyCriterion(convert(eltype(grid), -1e-4 * g_Earth)))
	validate_criterion_model(criterion, args...)
	MLD = MixedLayerDepth(criterion)
	return KernelFunctionOperation{Center, Center, Nothing}(MLD, grid, args...)
	end

Which actually returns a KernelFunctionOperation. Am I correct?

If so, that feels a bit confusing/convoluted to me. Is there a way to make this more straightforward?

@glwagner I think this overlaps with our discussion in #194, so maybe some points are still not connecting in my head...

[tomchor]

@jagoosw I addressed both mine and @glwagner's comments directly. Let me know if you don't agree with anything and I'll revert the commits.