Introduction

The coupling functionality in FABM was designed to enable local interaction between multiple model instances. Here, each instance is active in the same domain: in the pelagic, at the water surface, and/or at the bottom. However, numerous applications need the ability to couple depth-integrated variables (2D) with depth-explicit [pelagic] ones (3D). Notably, Eulerian higher trophic level (HTL) models such as ECOSMO E2E, SEAPODYM-LMTL and a spatially explicit implementation of the Community Size Spectrum Model represent HTL biomass by depth-integrated quantities that interact with an arbitrary part of the water column, for example, to obtain prey and dissolved oxygen, or to inject waste products such as CO₂. To underpin this interaction, the predator has a vertical “distribution”, “habitat” or “search range”, specified by weights that indicate what proportion of time the HTL spends in each layer. Such logic previously was difficult to implement in FABM: while possible, it required use of poorly documented APIs (e.g., do_column) and, due to the need for explicit and complex 2D-3D mapping, it led to hard-to-read code.

Version 2.1 of FABM added functionality to facilitate the mapping between the 3D pelagic environment (e.g., temperature, net primary production, prey) and 2D higher trophic levels, while respecting the vertical distribution of the predator and guaranteeing mass conservation. This is specifically designed to facilitate the coupling between lower and higher trophic level models.

This section describes the new functionality and also illustrates what happens behind the scenes. It is not essential to know the latter in order to work on a depth-integrated/higher-trophic level model. To get started quickly, take a look at the example included with FABM.

Base type

The primary way to use this new functionality is to have your predator module inherit from the new type_depth_integrated_particle type introduced in new Fortran module fabm_builtin_depth_mapping:

type, extends(type_depth_integrated_particle) :: type_depth_integrated_predator
   ...
end type

type_depth_integrated_particle provides subroutines to access to depth integrals and depth averages of depth-explicit dependencies. These routines can be used, for example, to access temperature or prey averaged over a predator's vertical habitat. The built-in depth integration logic only needs to be provided with weights for every layer, which define the depth distribution of your predator: to have the predator inhabit only the top 100 m of the water column, you would set these weights to 1 in cells at depths ≤ 100 m, and to 0 in all deeper cells. The weights are allowed to vary in time and space, and thus can account for behaviours such as diurnal or seasonal vertical migration. The weights are not restricted to either 0 or 1: they represent the relative fraction of time spent at each depth and can therefore take any real value. Formally, this fraction equals $\frac{w_i h_i}{\sum_k w_k h_k}$, with $w_i$ being the weight in layer $i$, and $h_i$ the layer's thickness. This expression also shows that only the relative value of the weights matters: scaling all weights $w$ with an arbitrary depth-independent constant has no impact on the fraction of time spent in a layer. Thus, to have the predator distribution track prey availability, the weights can simply be set equal to the prey concentration.

In many cases, you will want to create your own simple child model to calculate these weights as a function of environmental variables, depth and/or time. That allows you to prescribe custom recipes, e.g., “make the predator distribution proportional to pelagic prey but avoid areas where oxygen is below 5 mmol m-3”.

Dependencies on depth-averaged variables

To add a dependency on depth-averaged temperature to your predator, you add a surface (i.e., depth-independent) dependency to your type as usual:

type (type_surface_dependency_id) :: id_temp

This is registered in your initialize routine as usual:

call self%register_dependency(self%id_temp, 'temp', 'degrees_Celsius', 'depth-averaged temperature')

The one new step is to couple the new dependency to the actual depth-average of temperature, which is done with a single call to request_mapped_coupling:

call self%request_mapped_coupling(self%id_temp, standard_variables%temperature, average=.true.)

This sets up a child model that performs the depth-averaging, as shown in Fig 6. request_mapped_coupling will ensure that this child model will respect the vertical weights that specify the predator habitat. The argument average=.true. ensures that id_temp will receive the depth average $\frac{\sum_k w_k c_k h_k}{\sum_k w_k h_k}$, with $w_k$, $c_k$ and $h_k$ representing local distribution weights, temperature, and layer thickness, respectively. If the average argument is omitted or set to .false., id_temp receives the depth integral $\sum_k w_k c_k h_k$ instead.

Figure 6. A depth-integrated predator type with a dependency on depth-averaged temperature. All ellipses indicate variables; dependencies are indicated by ellipses with dashed border. The depth-integrated predator would have its own entry in fabm.yaml; conversely, the child model “integrator” is automatically created by the call to register_mapped_coupling.

The very similar request_mapped_coupling_to_model provides access to the depth average or depth integral of a variable from another model instance. For example, the following links id_prey_c to the depth-averaged total carbon from a coupled model instance named “prey”.

call self%request_mapped_coupling_to_model(self%id_prey_c, 'prey', standard_variables%total_carbon, average=.true.)

This “prey” instance must be coupled to an actual model instance in the coupling section of the predator in fabm.yaml. An example where the predator is coupled to a fixed-stoichiometry zooplankton instance “Z4” is shown in Fig 7.

Figure 7. A depth-integrated predator type with a dependency on depth-averaged prey carbon. All ellipses indicate variables; dependencies are indicated by ellipses with dashed border. Both the depth-integrated predator and the mesozooplankton would have its own entry in fabm.yaml; conversely, the child model “integrator” is automatically created by the call to register_mapped_coupling_to_model.

More information about coupling to model instances and the use of standard variables such as total_carbon can be found in Particles: grouped couplings and generalized access.

Distributing depth-integrated source terms over the pelagic

The above examples set up read-only dependencies: they allow you to easily obtain the (weighted) depth average or integral of a depth-explicit variable. What if you also want to return source terms for this variable? In that case, the syntax is the same, except that the receiving variable should be a state variable dependency (access: read, increment sources), rather than a plain dependency (read access only). For example, to link to a pelagic waste pool to which you want to direct a depth-integrated carbon flux, the dependency would be declared like:

type (type_surface_state_variable_id) :: id_waste_c

It would be registered in your initialize subroutine as usual:

call self%register_state_dependency(self%id_waste_c, 'waste_c', 'mmol C m-2', 'depth-integrated carbon waste')

Finally, it would be coupled to the depth-integrated total carbon of a “waste” instance with

call self%request_mapped_coupling_to_model(self%id_waste_c, 'waste', standard_variables%total_carbon)

In this case, we do not specify argument average=.true., since we will provide a depth-integrated waste source. Accordingly, the dependency is expressed per unit surface area, not per volume.

The above syntax is all nearly identical to the treatment of depth-averaged temperature and prey described earlier. However, the fact that id_waste_c is a state variable changes the behaviour of the depth-integrating child model, as shown in Fig 8. It now sets up a child instance that collects source terms for the depth-integrated carbon waste and redistributes them over the original depth-explicit waste; all while respecting the predator’s depth distribution.

Figure 8. A depth-integrated predator type with a dependency on depth-integrated carbon in a pelagic waste pool. Unlike previous examples (Figs 6, 7), this state variable dependency allows the predator to feed back a depth-integrated source term, which then gets redistributed over the original depth-explicit waste pool. All ellipses indicate variables; dependencies are indicated by ellipses with dashed border. Both the depth-integrated predator and the particulate matter would have its own entry in fabm.yaml; the child models “integrator” and “rate_distributor” are automatically created by the call to register_mapped_coupling_to_model.

The default rule for distributing a depth-integrated flux over a pelagic state variable is to direct fraction $\frac{w_i h_i}{\sum_k w_k h_k}$ of the depth-integrated flux $f$ to layer $i$. Here, $w$ represents the weights that define the predator's vertical distribution and $h$ the layer thickness. Expression, $\sum_k w_k h_k$ is the depth-integral of weights for the local water column. The local change in the pelagic state variable is the layer-integrated flux divided by layer thickness, i.e., $f \frac{w_i}{\sum_k w_k h_k}$. Note that this ensures that the flux is zero where the predator is absent (i.e., where $w_i=0$). Additionally, only the relative value of weights $w$ matter for the redistribution rule: scaling the weights with an arbitrary depth-independent constant does not change the final redistributed source term.

The rate distributor module supports one alternative redistribution rule, under which the flux directed into each layer is proportional to the local concentration of the pelagic state variable (as well as being proportional to predator weights $w$). Such proportionality between the local flux and the local concentration is appropriate for certain loss terms, notably the loss of prey due to predation. The local redirected flux is now proportional to $w_i c_i$, with $c_i$ denoting the local pelagic concentration. To ensure mass conservation (the depth integral of redirected flux must equal the originally specified depth-integrated flux), it follows that the local change must equal $f \frac{w_i c_i}{\sum_k w_k c_i h_k}$. To use this alternative redistribution rule, register_mapped_coupling_to_model must be called with argument proportional_change=.true.

Applying loss terms to pelagic prey

The final new component is the ability to access (loop over) all depth-integrated state variables of a coupled pelagic instance, and to provide depth-integrated rates of change for each. This is typically used to apply the same relative loss rate to all constituents of a prey (e.g., carbon, nitrogen, phosphorus, chlorophyll, etc.), without having to know the prey’s implementation details. It allows the same prey treatment as described near the end of section Particles: grouped couplings and generalized access. To use this functionality, you first add a model dependency to the model type:

type (type_model_id) :: id_prey_int

This is registered from initialize by calling the new register_mapped_model_dependency:

call self%register_mapped_model_dependency(self%id_prey_int, 'prey', proportional_change=.true., domain=domain_surface)

Here, the use of proportional_change=.true. specifies that when redistributing prey loss terms over the vertical, the local change in prey must be proportional to the local prey concentration (see the discussion of the alternative distribution rule above). By specifying domain=domain_surface, we ensure that all depth-integrated prey variables become available as id_prey_int%surface_state. These variables can then be read and changed by specifying source terms from do_surface:

do istate = 1, size(self%id_prey_int%surface_state)
   _GET_SURFACE_(self%id_prey_int%surface_state(istate), p)
   _ADD_SURFACE_SOURCE_(self%id_prey_int%surface_state(istate), -prey_loss_rate * p)
end do

Here, prey_loss_rate is the specific loss rate of the depth-integrated prey, in s-1. In the background, registered_mapped_model_dependency takes care of setting up depth-integrated variables for each of the prey's original pelagic state variables, along with the appropriate integrator child models and rate distributors, as shown in Fig. 9. It is worth emphasizing, however, that all new functionality can be used without being aware of the underlying implementation.

Figure 9. A depth-integrated predator type with a dependency on a pelagic prey instance. By calling register_mapped_model_dependency, the predator gets access to all depth-integrated state variables of the prey, irrespective of their identity. It is then free to prescribe rates of change for each, e.g., to apply the same specific loss rate. These rates get redistributed over the original depth-explicit state variables of the prey. All ellipses indicate variables; dependencies are indicated by ellipses with dashed border. Both the depth-integrated predator and the mesozooplankton would have its own entry in fabm.yaml; the child models “particle_integrator” and all necessary instances of “integrator” and “rate_distributor” are automatically created by the call to register_mapped_model_dependency.

A comprehensive example that uses combines all functionality described above is included in FABM.