Implement surface aerodynamcis following standard Monin-Obukhov similarity theory

[bgroenks96]

Currently, the SurfaceEnergyBalance relies on a simple parameterization ConstantAerodynamics that provides only a constant value for the surface drag coefficient. This is quite far from realistic since, in reality, this coefficient varies considerably based on transient atmospheric conditions and surface properties.

I would suggest that we start by migrating the implementation of Monin-Obukhov in CryoGrid.jl or by implementing the scheme from PALADYN/CLIMBER-X. Both would probably be good enough for now.

@olivierbonte Do you want to maybe take a crack at this?

CC @milankl

[bgroenks96]

Also @glwagner maybe you want to weigh in on whether or not this falls within the scope of NumericalEarth.jl. If NumericalEarth becomes the main engine for coupling models, then we could conceivably just have EarthSystemModel(atmosphere = SpeedyAtmosphere(...), land = TerrariumLand(...)) and then NumerialEarth.jl would host the more complex boundary layer aerodynamics?

Edit: I guess this doc page basically answers my question.

[glwagner]

You can use NumericalEarth to compute fluxes between a Terrarium land model and a prescribed atmosphere (for example). This is how we do ocean simulations with MO flux models.

The proposal of NumericalEarth is to put flux models in one place (EarthSystemModel) rather than duplicate their implementation in different components.

To illustrate why this is powerful we are building a suite of examples with a prognostic atmosphere over various surfaces:

https://numericalearth.github.io/NumericalEarthDocumentation/dev/literated/breeze_over_two_oceans/

right now we have a slab ocean and hydrostatic ocean, but a future PR will add prescribed SST and nonhydrostatic ocean.

Is there something missing from the NumericalEarth solvers that you need for Terrarium? You could consider building features for the EarthSystemModel flux interface.

[olivierbonte]

I'd be interested to give it a shot! I see that the docs of v0.1 are taking shape (#84), so once I go over this, I am hopefully in an good spot to give this a shot :) I'll probably have time for this next week.

Having a quick look at the NumericalEarth docs page on this, it seems like a lot of useful stuff is indeed already in place, although it would probably need some extensions for land specific stuff e.g. surrounding roughness lengths and displacements height for vegetated land surfaces.

I am also not that familiar with the formulation outlined in https://numericalearth.github.io/NumericalEarthDocumentation/stable/interface_fluxes/#Buoyancy-flux-and-stability-of-the-near-surface-atmosphere, most implementations I have seen rely on the Bulk Richardson number (see e.g. Chapter 5 of the CLM documentation).

For non-vegetated areas, the soil temperature of the first layer can be used to calculate the delta in temperature (see CLM5 docs). Vegetated surfaces bring extra complexity of having unknown temperatures of both the canopy surface and the canopy air (see CLM5 docs again).

So I am definitely interested to help, but I think I'd need some extra guidance on where you'd want the implementation (Terrarium.jl and/or NumericalEarth.jl), what exact formulation you are interested in (CyroGrid, Paladyn, CLM...) and how extensive you want this to be (bare soil only, including snow, including vegetation...)

[milankl]

Btw @gregrmunday is developing a NN-based surface roughness scheme in SpeedyWeather/SpeedyWeather.jl#959 which seems to pass the Pangaea stress test too. Maybe helpful!

[bgroenks96]

So I am definitely interested to help, but I think I'd need some extra guidance on where you'd want the implementation (Terrarium.jl and/or NumericalEarth.jl), what exact formulation you are interested in (CyroGrid, Paladyn, CLM...) and how extensive you want this to be (bare soil only, including snow, including vegetation...)

I am a bit unsure on this myself. This largely hinges on whether or not we want Terrarium.jl to still be somewhat standalone and enable, for example, offline simulations with prescribed meteorological inputs. I lean towards saying yes we do want that because it seems somehow unreasonable to expect the user to install NumericalEarth and set up an EarthSystemModel just to do pure land simulations. I know @glwagner sees this differently.

On the other hand, there shouldn't be major differences between the surface flux calculations for a prescribed vs. prognostic atmosphere. So this strategy would inevitably lead to some code duplication; there probably already is some since Terrarium currently has a basic SEB already. This would be an argument in favor of relying on NumericalEarth for all flux coupling regardless of the offline vs. online case.

My main concern @glwagner is that a lot of surface flux calculations are very tightly coupled with land surface properties. The resistance terms directly depend on skin temperature and land parameters, and if you do this the "right" way, the full vegetation energy module needs to be re-evaluated at each iteration of the solver for skin temperature. Similarly, as we discussed before, land models do not interface with the atmosphere only via fluxes; vegetation models directly consume other atmospheric state variables like air temperature and humidity to compute stress factors. So I think we wold basically end up needing to move a lot of the land physics (including entire chunks of the vegetation module) to NumericalEarth.jl and then it would be unclear what the role of Terrarium would even be as a standalone package.

[olivierbonte]

Agree with what you're saying here @bgroenks96. As a user, I think being able to do offline runs e.g. at an eddy covariance site using only Terrarium.jl would be easiest.

Conceptually, it is indeed the cleanest (and most correct) to consider variables like surface temperature or temperature in the canopy air space as things to solve for in the nonlinear solve that also includes the calculation of the surface fluxes. In SUMMA (modular hydrology framework, https://doi.org/10.1029/2024MS004256), this is also why they moved to a differential algebraic equation solver (SUNDIALS) to solve both the ODE and the nonlinear algebraic constraints (mainly related to the surface energy balance) together.

[glwagner]

My main concern @glwagner is that a lot of surface flux calculations are very tightly coupled with land surface properties. The resistance terms directly depend on skin temperature and land parameters, and if you do this the "right" way, the full vegetation energy module needs to be re-evaluated at each iteration of the solver for skin temperature. Similarly, as we discussed before, land models do not interface with the atmosphere only via fluxes; vegetation models directly consume other atmospheric state variables like air temperature and humidity to compute stress factors. So I think we wold basically end up needing to move a lot of the land physics (including entire chunks of the vegetation module) to NumericalEarth.jl and then it would be unclear what the role of Terrarium would even be as a standalone package.

I think this is a good argument to implement surface flux calculations only once. Otherwise, you have to implement these tightly coupled formulations in two places.

As a thought experiment --- say that I want to develop a new flux parameterization for atmosphere-land simulations. Should I implement this in Terrarium or NumericalEarth? If I am interested in prognostic atmospheres (this is my case), then I probably will use NumericalEarth. But then this parameterization is not available in Terrarium without porting it. Conversely, if I am interested in land-only simulations, then perhaps I should use Terrarium. But then, it will not be available in NumericalEarth without porting. In the worst case, it doubles the work. Also, the user interface for Terrarium and NumericalEarth can drift.

the full vegetation energy module needs to be re-evaluated at each iteration of the solver for skin temperature.
Similarly, as we discussed before, land models do not interface with the atmosphere only via fluxes; vegetation models directly consume other atmospheric state variables like air temperature and humidity to compute stress factors.

If I understand correctly, this algorithm must be supported in NumericalEarth for coupled atmosphere-land simulations. Otherwise, we cannot perform coupled atmosphere-land simulations with NumericalEarth.

So I think we wold basically end up needing to move a lot of the land physics (including entire chunks of the vegetation module) to NumericalEarth.jl

The task is not to move things but to have interfaces for them in NumericalEarth. Terrarium will then use those interfaces (it will be implemented in NumericalEarthTerrariumExt, per the design we currently use for supporting different components).

and then it would be unclear what the role of Terrarium would even be as a standalone package.

In the case that NumericalEarth supports more than one land model (eg ClimaLand and Terrarium), the role of Terrarium as an independent package is quite clear I think.

Another possibility is to implement a new package for surface flux calculations and use this in both NumericalEarth and Terrarium. However, given how the surface flux implementation is tightly coupled to the land model itself, I feel this isn't very different in practice from simply implementing it in NumericalEarth. However it is the approach that SurfaceFluxes.jl takes.

I lean towards saying yes we do want that because it seems somehow unreasonable to expect the user to install NumericalEarth and set up an EarthSystemModel just to do pure land simulations

The design in NumericalEarth is to build bespoke aliases and constructors for specific cases. For example, with a land Simulation and PrescribedAtmosphere, we have a type alias like LandModel or LandOnlyModel.

I am not sure it is possible to evaluate if it is "somehow unreasonable" without having example code. Eg -- what specifically is unreasonable? If the user interface resembles what we do for the ocean model then it might look like

using NumericalEarth
using Terrarium

grid = RectilinearGrid(size=16, x=longitude, y=latitude, z=(-10, 0), topology=(Flat, Flat, Bounded))
land = land_simulation(grid)
atmosphere = ERA5PrescribedAtmosphere(grid)
model = LandOnlyModel(land; atmosphere)

we should look at the corresponding user interface with your Terrarium-only design.

Also, if you have two different implementations for surface fluxes, I think it is very importantt that the flux parameterizations can be configured identically in the two codebases.

[bgroenks96]

As a thought experiment --- say that I want to develop a new flux parameterization for atmosphere-land simulations. Should I implement this in Terrarium or NumericalEarth? If I am interested in prognostic atmospheres (this is my case), then I probably will use NumericalEarth. But then this parameterization is not available in Terrarium without porting it. Conversely, if I am interested in land-only simulations, then perhaps I should use Terrarium. But then, it will not be available in NumericalEarth without porting. In the worst case, it doubles the work. Also, the user interface for Terrarium and NumericalEarth can drift.

This seems to presuppose that the surface flux computations are always the same for prognostic and prescribed atmospheres. Is that the case? I really don't know because most of my experience is in implementing surface fluxes offline (prescribed forcing data).

If I understand correctly, this algorithm must be supported in NumericalEarth for coupled atmosphere-land simulations. Otherwise, we cannot perform coupled atmosphere-land simulations with NumericalEarth.

There can and maybe should be different schemes available with different levels of complexity. A full implicit solve of the skin temperature might be nice to have in a land-only simulation but might be too expensive and not worth it for a coupled simulation. PALADYN avoids the iterative solve by doing a single iteration based on a Taylor expansion, for example.

The task is not to move things but to have interfaces for them in NumericalEarth. Terrarium will then use those interfaces (it will be implemented in NumericalEarthTerrariumExt, per the design we currently use for supporting different components).

But how exactly will this work? Terrarium is not going to depend on NumericalEarth, right? So if the only way to actually do a simulation including surface fluxes is with NumericalEarth, this makes Terrarium essentially non-functional as a standalone package, at least for full land model configurations.

Maybe that's ok and we can just provide minimally functional approximations based on common simplifications like Penman-Monteith, but it's something worth thinking about whether or not that's really mutually beneficial.

In the case that NumericalEarth supports more than one land model (eg ClimaLand and Terrarium), the role of Terrarium as an independent package is quite clear I think.

Not really. ClimaLand is a fully functional standalone land model and they do currently duplicate the functionality we're discussing. As I understand it, you're arguing against further duplication, which I sympathize with, but you can't really compare Terrarium to ClimaLand in this case.

Another possibility is to implement a new package for surface flux calculations and use this in both NumericalEarth and Terrarium. However, given how the surface flux implementation is tightly coupled to the land model itself, I feel this isn't very different in practice from simply implementing it in NumericalEarth. However it is the approach that SurfaceFluxes.jl takes.

Is there any reason why we can't just use SurfaceFluxes.jl? It looks much better now than the last time I looked at it ~6 months ago.

we should look at the corresponding user interface with your Terrarium-only design.

I think your snippet is a bit misleading. Clearly land_simulation cannot be initialized with only a grid... the user has to provide a Terrarium model or simulation at the very least. But something comparable would be our ERA5-Land example:

# Load land-sea mask at ~1° resolution
land_sea_frac = convert.(Float32, dropdims(Raster("inputs/era5-land_land_sea_mask_N72.nc"), dims = Ti))
land_sea_frac_field = RingGrids.FullGaussianField(Matrix(land_sea_frac), input_as = Matrix)
heatmap(land_sea_frac_field)


# Load ERA5-Land 2 meter air temperature at ~1° resolution
Tair_raster = Raster("inputs/external/era5-land/2m_temperature/era5_land_2m_temperature_2023_N72.nc")

# Set up grids
land_mask = land_sea_frac_field .> 0.5 # select only grid points with > 50% land
Nz = 30 # number of soil layers
grid = ColumnRingGrid(arch, ExponentialSpacing(N = Nz), land_mask)

# Construct input sources
Tair_forcing = InputSource(grid, rebuild(Tair_raster, name = :Tair))
model = SoilModel(grid)
boundary_conditions = PrescribedSurfaceTemperature(:Tair)
Tsurf_0 = Tair_raster[Ti(1)][findall(land_mask)]
# Initial conditions
initializers = (
    # steady-ish state initial condition for temperature
    temperature = (x, z) -> Tsurf_0[Int(round(x)) + 1] - 0.02 * z,
    # dry soil
    saturation_water_ice = 1.0,
)
integrator = initialize(model, ForwardEuler(), Tair_forcing; initializers, boundary_conditions)

This is more code because it doesn't rely so much on defaults... and the interface of course still needs further refinement; e.g. we don't currently have a specialized implementation of AbstractAtmsophere for ERA5 (though I am not sure I understand why that is necessary?).

Also, if you have two different implementations for surface fluxes, I think it is very importantt that the flux parameterizations can be configured identically in the two codebases.

Yes, I agree that is very important.

[glwagner]

So if the only way to actually do a simulation including surface fluxes is with NumericalEarth, this makes Terrarium essentially non-functional as a standalone package, at least for full land model configurations.

It's the same for Oceananigans and Breeze. They cannot be used for OMIP or AMIP like configurations. They can only be used for idealized cases with either prescribed fluxes or simple flux models, like constant coefficients.

This seems to presuppose that the surface flux computations are always the same for prognostic and prescribed atmospheres. Is that the case? I really don't know because most of my experience is in implementing surface fluxes offline (prescribed forcing data).

I don't see why this matters. NumericalEarth supports both prescribed and prognostic atmospheres.

There can and maybe should be different schemes available with different levels of complexity. A full implicit solve of the skin temperature might be nice to have in a land-only simulation but might be too expensive and not worth it for a coupled simulation.

We use a full implicit skin temperature solve for sea ice now and the plan is to also use this for ocean as well. There are lots of paths towards optimization (tabulating expensive functions, low precision, etc). I think we will be able to afford it.

But how exactly will this work? Terrarium is not going to depend on NumericalEarth, right?

We will design it based on the requirements. Sea ice and snow have this requirement: to compute the surface energy balance, we need to be able to compute the "interior" heat flux that is coming from the bulk to the surface. The surface temperature is then determined by a balance between all of the fluxes. The surface temperature solver can call a function to retrieve the interior heat flux, which can invoke the land model as needed. We can also support a variety or hierarchy of interfaces. For example, if the heat flux is based on a diffusivity, a bespoke interface for diffusive heat fluxes can be provided (in which the two inputs are diffusivity and bulk temperature).

This design is crucial: it creates a situation in which NumericalEarth's domain of interest includes surface AND near-surface processes. For example, to compute atmosphere-land fluxes, NumericalEarth's design must be compatible with the latest concepts for modeling near-surface land fluxes. Hopefully this will help eliminate inconsistent modeling of interfaces which is a problem with some legacy coupled models.

[glwagner]

Not really. ClimaLand is a fully functional standalone land model and they do currently duplicate the functionality we're discussing. As I understand it, you're arguing against further duplication, which I sympathize with, but you can't really compare Terrarium to ClimaLand in this case.

I am not comparing Terrarium to ClimaLand. I said "In the case that NumericalEarth supports more than one land model (eg ClimaLand and Terrarium), the role of Terrarium as an independent package is quite clear I think." This has nothing to do with comparing Terrarium and ClimaLand.

[glwagner]

Is there any reason why we can't just use SurfaceFluxes.jl? It looks much better now than the last time I looked at it ~6 months ago.

No, there is no reason why you cannot use it.

[glwagner]

I think your snippet is a bit misleading. Clearly land_simulation cannot be initialized with only a grid... the user has to provide a Terrarium model or simulation at the very least.

In this design land_simulation builds a Terrarium simulation. We also have ocean_simulation(grid). Why is this misleading? It is actually the interface we use:

https://github.com/NumericalEarth/NumericalEarth.jl/blob/558bf7605a3ef50277f32ea8bd1d0e89c0fabdf2/examples/near_global_ocean_simulation.jl#L78

This is more code because it doesn't rely so much on defaults...

All of that complexity is orthogonal to whether you use LandOnlyModel or not. Using NumericalEarth adds at most one additional line to build LandOnlyModel. This is the point.

we don't currently have a specialized implementation of AbstractAtmsophere for ERA5 (though I am not sure I understand why that is necessary?).

I hope that you do not need to have such an abstraction, and that you can simply use the representations in NumericalEarth.