# Combined physical and biogeochemical assessment of mesoscale eddy parameterisations in ocean models: Eddy-induced advection at eddy-permitting resolutions

X Ruan, Damien Couespel, Marina Lévy, J Li, J Mak, Y Wang

► To cite this version:

X Ruan, Damien Couespel, Marina Lévy, J Li, J Mak, et al.. Combined physical and biogeochemical assessment of mesoscale eddy parameterisations in ocean models: Eddy-induced advection at eddy-permitting resolutions. Ocean Modelling, 2024, 190, pp.102396. 10.1016/j.ocemod.2024.102396. hal-04620284

HAL Id: hal-04620284

https://hal.science/hal-04620284

Submitted on 21 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Combined physical and biogeochemical assessment of mesoscale eddy parameterisations in ocean models: Eddy-induced advection at eddy-permitting resolutions

Xi Ruan, Damien Couespel, Marina Lévy, Jiying Li, Julian Mak, Yan Wang

► To cite this version:

Xi Ruan, Damien Couespel, Marina Lévy, Jiying Li, Julian Mak, et al.. Combined physical and biogeochemical assessment of mesoscale eddy parameterisations in ocean models: Eddy-induced advection at eddy-permitting resolutions. Ocean Modelling, 2024, 190, pp.102396. 10.1016/j.ocemod.2024.102396. hal-04620284v2

HAL Id: hal-04620284

https://hal.science/hal-04620284v2

Submitted on 21 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

# Combined physical and biogeochemical assessment of mesoscale eddy parameterisations in ocean models: Eddy-induced advection at eddy-permitting resolutions

X. Ruan<sup>a,*</sup>, D. Couespel<sup>b,c</sup>, M. Lévy<sup>c</sup>, J. Li<sup>a,d</sup>, J. Mak<sup>a,e,**</sup>, Y. Wang<sup>a,d</sup>

## ARTICLE INFO

Dataset link: http://dx.doi.org/10.5281/zenod
0.11498192

Keywords:

Scale-aware mesoscale eddy parameterisation

Ocean circulation

Climate change

Numerical modelling

Biophysical interactions

## ABSTRACT

Ocean general circulation models at the eddy-permitting regime are known to under-resolve the mesoscale eddy activity and associated eddy-mean interaction. Under-resolving the mesoscale eddy field has consequences for the resulting mean state, affecting the modelled ocean circulation and biogeochemical responses, and impacting the quality of climate projections. There is an ongoing debate on whether and how a parameterisation should be utilised in the eddy-permitting regime. Focusing on the Gent-McWilliams (GM) based parameterisations, it is known that, on the one hand, not utilising a parameterisation leads to insufficient eddy feedback and results in biases. On the other hand, utilising a parameterisation leads to double-counting of the eddy feedback, and introduces other biases. A recently proposed approach, known as splitting, modifies the way GM-based schemes are applied in eddy-permitting regimes, and has been demonstrated to be effective in an idealised Southern Ocean channel model. In this work, we evaluate whether the splitting approach can lead to improvements in the physical and biogeochemical responses in an idealised double gyre model. Compared with a high resolution mesoscale eddy resolving model truth, the use of the GM-based GEOMETRIC parameterisation together with splitting in the eddy-permitting regime leads to broad improvements in the control pre-industrial scenario and an idealised climate change scenario, over models with and models without the GM-based GEOMETRIC parameterisation active. While there are still some deficiencies, particularly in the subtropical region where the transport is too weak and may need momentum re-injection to reduce the biases, the present work provides further evidence in support of using the splitting procedure together with a GM-based parameterisation in ocean general circulation models at eddy-permitting resolutions.

## 1. Introduction

Mesoscale eddy effects are essential for shaping the ocean circulation, marine ecosystems and global climate via the associated eddy-mean interaction (e.g., Griffies et al., 2015; Fox-Kemper et al., 2019; Beech et al., 2022). One such effect is the slumping of isopycnals that is normally associated with baroclinic instability, which releases the large-scale available potential energy (from input via Ekman processes, buoyancy forcing or otherwise) into small-scale eddy energy (e.g., Gent and McWilliams, 1990; Gent et al., 1995), with

associated impacts on the mean stratification. Another is the additional diffusion of tracers, such as thermodynamic or biogeochemical variables, along isoneutral directions (e.g., Redi, 1982; Griffies, 1998; Jones and Abernathey, 2019; Holmes et al., 2022), which can affect the ventilation rate of tracers (although there are recent works that suggest such isoneutral diffusion can also have a non-negligible effect on shaping the stratification; e.g., Chouksey et al. 2022). Yet another is the re-injection of eddy energy back into larger-scales and forcing the larger-scale circulation (e.g., Charney, 1971; Waterman and Jayne,

<footnote><sup>a</sup> Department of Ocean Science, Hong Kong University of Science and Technology, Hong Kong</footnote>

<footnote><sup>b</sup> NORCE Norwegian Research Centre, Bjerknes Centre for Climate Research, Norway</footnote>

<footnote><sup>c</sup> Sorbonne Université, LOCEAN-IPSL, CNRS/IRD/MNHN, France</footnote>

<footnote><sup>d</sup> Center for Ocean Research in Hong Kong and Macau, Hong Kong University of Science and Technology, Hong Kong</footnote>

<footnote><sup>e</sup> National Centre of Oceanography, Southampton, United Kingdom</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>** Corresponding author at: Department of Ocean Science, Hong Kong University of Science and Technology, Hong Kong.
E-mail addresses: xruanaa@connect.ust.hk (X. Ruan), julian.c.l.mak@googlemail.com (J. Mak).</footnote>
<figure><img src="image_2.png" /></figure>

# Combined physical and biogeochemical assessment of mesoscale eddy parameterisations in ocean models: Eddy-induced advection at eddy-permitting resolutions

X. Ruan<sup>a,*</sup>, D. Couespel<sup>b,c</sup>, M. Lévy<sup>c</sup>, J. Li<sup>a,d</sup>, J. Mak<sup>a,e,**</sup>, Y. Wang<sup>a,d</sup>

## ARTICLE INFO

Dataset link: http://dx.doi.org/10.5281/zenod
0.11498192

Keywords:

Scale-aware mesoscale eddy parameterisation

Ocean circulation

Climate change

Numerical modelling

Biophysical interactions

## ABSTRACT

Ocean general circulation models at the eddy-permitting regime are known to under-resolve the mesoscale eddy activity and associated eddy-mean interaction. Under-resolving the mesoscale eddy field has consequences for the resulting mean state, affecting the modelled ocean circulation and biogeochemical responses, and impacting the quality of climate projections. There is an ongoing debate on whether and how a parameterisation should be utilised in the eddy-permitting regime. Focusing on the Gent-McWilliams (GM) based parameterisations, it is known that, on the one hand, not utilising a parameterisation leads to insufficient eddy feedback and results in biases. On the other hand, utilising a parameterisation leads to double-counting of the eddy feedback, and introduces other biases. A recently proposed approach, known as splitting, modifies the way GM-based schemes are applied in eddy-permitting regimes, and has been demonstrated to be effective in an idealised Southern Ocean channel model. In this work, we evaluate whether the splitting approach can lead to improvements in the physical and biogeochemical responses in an idealised double gyre model. Compared with a high resolution mesoscale eddy resolving model truth, the use of the GM-based GEOMETRIC parameterisation together with splitting in the eddy-permitting regime leads to broad improvements in the control pre-industrial scenario and an idealised climate change scenario, over models with and models without the GM-based GEOMETRIC parameterisation active. While there are still some deficiencies, particularly in the subtropical region where the transport is too weak and may need momentum re-injection to reduce the biases, the present work provides further evidence in support of using the splitting procedure together with a GM-based parameterisation in ocean general circulation models at eddy-permitting resolutions.

## 1. Introduction

Mesoscale eddy effects are essential for shaping the ocean circulation, marine ecosystems and global climate via the associated eddy-mean interaction (e.g., Griffies et al., 2015; Fox-Kemper et al., 2019; Beech et al., 2022). One such effect is the slumping of isopycnals that is normally associated with baroclinic instability, which releases the large-scale available potential energy (from input via Ekman processes, buoyancy forcing or otherwise) into small-scale eddy energy (e.g., Gent and McWilliams, 1990; Gent et al., 1995), with

associated impacts on the mean stratification. Another is the additional diffusion of tracers, such as thermodynamic or biogeochemical variables, along isoneutral directions (e.g., Redi, 1982; Griffies, 1998; Jones and Abernathey, 2019; Holmes et al., 2022), which can affect the ventilation rate of tracers (although there are recent works that suggest such isoneutral diffusion can also have a non-negligible effect on shaping the stratification; e.g., Chouksey et al. 2022). Yet another is the re-injection of eddy energy back into larger-scales and forcing the larger-scale circulation (e.g., Charney, 1971; Waterman and Jayne,

<footnote><sup>a</sup> Department of Ocean Science, Hong Kong University of Science and Technology, Hong Kong</footnote>

<footnote><sup>b</sup> NORCE Norwegian Research Centre, Bjerknes Centre for Climate Research, Norway</footnote>

<footnote><sup>c</sup> Sorbonne Université, LOCEAN-IPSL, CNRS/IRD/MNHN, France</footnote>

<footnote><sup>d</sup> Center for Ocean Research in Hong Kong and Macau, Hong Kong University of Science and Technology, Hong Kong</footnote>

<footnote><sup>e</sup> National Centre of Oceanography, Southampton, United Kingdom</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>** Corresponding author at: Department of Ocean Science, Hong Kong University of Science and Technology, Hong Kong.
E-mail addresses: xruanaa@connect.ust.hk (X. Ruan), julian.c.l.mak@googlemail.com (J. Mak).</footnote>
2012; Waterman and Hoskins, 2013; Bachman, 2019). The representation of such processes in numerical ocean models can have a leading order impact on the resulting model responses (e.g., Fox-Kemper et al., 2019). Thus, it is important to quantify the impacts of such eddy effects on the model responses, which constrains future model projections and their dependence on eddy processes, and helps inform management and policy decisions (e.g., Hoegh-Guldberg et al., 2014; Hewitt et al., 2020).

Most ocean models employed thus far for Earth System Models are coarse resolution models, in which the models do not permit an explicit representation of ocean mesoscale eddy effects. Eddy effects are often parameterised in these models, such as by the Gent-McWilliams (GM) based schemes to mimic the slumping of isopycnals (e.g., Gent and McWilliams, 1990; Gent et al., 1995; Eden and Greatbatch, 2008; Marshall et al., 2012; Mak et al., 2022b), diffusion of tracers along the isoneutral direction (e.g., Redi, 1982; Griffies, 1998; Holmes et al., 2022), momentum backscatter to re-inject eddy energy (e.g., Bachman, 2019; Jansen et al., 2019; Juricke et al., 2020; Yankovsky et al., 2023), or possibly related machine learning approaches (e.g., Perezhogin et al., 2023). While many works consider the physical or biogeochemical response to such parameterisations (e.g., Pradal and Gnanadesikan, 2014; Berthet et al., 2019; Swearer et al., 2019; Séférían et al., 2019; Mak et al., 2022b), of particular relevance to the present work are those of Couespel et al. (2021) and Ruan et al. (2023), who consider the joint physical and biogeochemical response to parameterisation choices. Those works, although within the context of idealised ocean models and simplified climate change scenarios, comprehensively assess both responses in coarse resolution models, highlighting how some choices of GM-based parameterisations are able to improve aspects of the model responses and sensitivities. In particular, the work of Ruan et al. (2023) highlights a case where one can obtain a reasonable biogeochemical response but for an essentially inconsistent physical response, i.e., where we get a 'right' response but not necessarily for the 'right' reasons.

The present work builds upon the previous work of Ruan et al. (2023) in coarse resolution models, by performing an analogous examination in the eddy-permitting regime. With increasing computational power, it is increasingly feasible to allow for an explicit representation of mesoscale eddies in numerical ocean models, even in global configuration Earth System Models. While it is known that increasing the spatial resolution of the ocean models can lead to a reduction of model biases (e.g., Roberts et al., 2020; Hewitt et al., 2020, 2022; Beech et al., 2022), other issues arise, particularly as to whether mesoscale eddy parameterisations should be dispensed with or employed, and if employed, how they are employed. Without parameterisations, it is known that the eddy processes are mis-represented, and the associated eddy-mean feedback is too weak. On the other hand, it is known that GM-based parameterisations utilised as-is end up damping explicit eddy fluctuations. One suggestion then is to switch off the GM-based parameterisations when the model is regarded as eddy-permitting, via the use of a resolution function (e.g., Hallberg, 2013). Another is to accept there is some damping by GM-based schemes, but re-energise via backscatter approaches (e.g., Jansen et al., 2015b,a, 2019; Bachman, 2019). A proposal considers anisotropic versions of the GM-based parameterisation (e.g., Smith and Gent, 2004). Some works advocate for backscatter only (e.g., Juricke et al., 2020; Chang et al., 2023; Yankovsky et al., 2023), arguing that backscatter re-energises the explicit eddy activity and catalyses for the extra eddy-mean interaction.

In the present work, the principal focus is on quantifying the benefits and/or deficiencies arising from a procedure suggested by Mak et al. (2023) termed field splitting, which modifies the way GM-based schemes are applied. No resolution function, anisotropic GM, backscatter or isoneutral diffusion is employed in this work, although the splitting procedure is not necessarily mutually exclusive of those modelling choices. We apply the splitting procedure to the model and experimental procedure of Ruan et al. (2023) in the eddy-permitting regime, where we expect improvements for reasons to be detailed.

The present work further tests the splitting procedure for a model in a different ocean-relevant setting (a representative mid-latitude gyre system), and comprehensively assesses both the modelled physical and biogeochemical responses, to inform future works using more realistic physical settings, biogeochemical models, and/or forcing scenarios of more direct relevance to climate projection exercises.

In Section 2, we provide the technical problem statement relating to the use of a GM-based parameterisation in an eddy-permitting regime, as well as the motivation and overview for the splitting procedure of Mak et al. (2023). We leverage the model and experimental procedure based largely on the previous work of Ruan et al. (2023) to test our scientific hypotheses, but with the caveat that the results we present in the main body of the work do not employ any isoneutral diffusion, for reasons to be elaborated upon. In Section 3 we provide a detailed comparison of the physical and biogeochemical responses in the set of models under a control pre-industrial scenario. In Section 4 we consider the analogous responses under an idealised climate change scenario, and additionally assess the associated sensitivities. The paper closes with conclusions and discussions in Section 5, detailing implications, limitations and outlooks in light of the present results. In Appendix, we reiterate and further elaborate on why the present work does not include isoneutral diffusion, and provide sample numerical evidence on why there are subtleties with the use of isoneutral diffusion together with state-aware GM-based parameterisations, such as that employed in this work.

## 2. Problem statement, field splitting approach, and numerical set up

### 2.1. Problem statement

The underlying problem in this work relates to whether and how a GM-based parameterisation should be used when the model allows for an explicit representation of mesoscale eddies. Consider a Reynolds decomposition of the velocity (specifically, the advective velocity in the tracer equation only)

$$u = \bar{u} + u' + u^*, \quad \bar{u}' = 0, \qquad (1)$$

where an overbar denotes a Reynolds average (time-average is considered in this work), a prime denotes a deviation from that average such that the average of the deviation is zero (so $u'$ is associated with the explicit eddies), and a star denotes any parameterised component we may wish to add on (so $u^*$ is associated with the parameterised eddies). In the eddy-rich/resolving case without a parameterisation, we would set $u^* = 0$. In a coarse resolution calculation, $u'$ would be effectively zero, and we might mimic the effect of the missing $u'$ by $u^*$, such as via the GM specification (e.g., Gent et al., 1995)

$$u^* = \nabla \times (e_z \times \kappa_{gm} s), \quad s = \frac{\nabla_H \rho}{-\partial \rho / \partial z}. \qquad (2)$$

where $\nabla$ is the gradient operator ($\nabla_H$ is the gradient operator only in the horizontal), $e_z$ is the unit vector pointing in the vertical, $\kappa_{gm}$ is the GM or eddy-induced velocity coefficient (in units of $m^2 s^{-1}$), $s$ is the vector encoding the isopycnal slopes in the horizontal directions, and $\rho$ denotes the relevant density variable of interest.

A modelling problem in the eddy-permitting regime is that $u'$, while weak, is not necessarily negligible, and the question is whether $u^*$ should be included or not. If $u^* = 0$ (by setting $\kappa_{gm} = 0$ for example), then the explicit eddies are too weak, leading to a rather weak feedback onto the mean state. As such, we might expect the resulting baroclinic mean flow to be too strong, associated for example with an overly deep stratification (e.g., Fig. 7c of Mak et al., 2023), resulting in an overly strong meridional overturning circulation. A larger lateral transport in the present double-gyre system to be investigated would be expected to lead to meridional heat transport with positive biases, while the
2012; Waterman and Hoskins, 2013; Bachman, 2019). The representation of such processes in numerical ocean models can have a leading order impact on the resulting model responses (e.g., Fox-Kemper et al., 2019). Thus, it is important to quantify the impacts of such eddy effects on the model responses, which constrains future model projections and their dependence on eddy processes, and helps inform management and policy decisions (e.g., Hoegh-Guldberg et al., 2014; Hewitt et al., 2020).

Most ocean models employed thus far for Earth System Models are coarse resolution models, in which the models do not permit an explicit representation of ocean mesoscale eddy effects. Eddy effects are often parameterised in these models, such as by the Gent-McWilliams (GM) based schemes to mimic the slumping of isopycnals (e.g., Gent and McWilliams, 1990; Gent et al., 1995; Eden and Greatbatch, 2008; Marshall et al., 2012; Mak et al., 2022b), diffusion of tracers along the isoneutral direction (e.g., Redi, 1982; Griffies, 1998; Holmes et al., 2022), momentum backscatter to re-inject eddy energy (e.g., Bachman, 2019; Jansen et al., 2019; Juricke et al., 2020; Yankovsky et al., 2023), or possibly related machine learning approaches (e.g., Perezhogin et al., 2023). While many works consider the physical or biogeochemical response to such parameterisations (e.g., Pradal and Gnanadesikan, 2014; Berthet et al., 2019; Swearer et al., 2019; Séférían et al., 2019; Mak et al., 2022b), of particular relevance to the present work are those of Couespel et al. (2021) and Ruan et al. (2023), who consider the joint physical and biogeochemical response to parameterisation choices. Those works, although within the context of idealised ocean models and simplified climate change scenarios, comprehensively assess both responses in coarse resolution models, highlighting how some choices of GM-based parameterisations are able to improve aspects of the model responses and sensitivities. In particular, the work of Ruan et al. (2023) highlights a case where one can obtain a reasonable biogeochemical response but for an essentially inconsistent physical response, i.e., where we get a 'right' response but not necessarily for the 'right' reasons.

The present work builds upon the previous work of Ruan et al. (2023) in coarse resolution models, by performing an analogous examination in the eddy-permitting regime. With increasing computational power, it is increasingly feasible to allow for an explicit representation of mesoscale eddies in numerical ocean models, even in global configuration Earth System Models. While it is known that increasing the spatial resolution of the ocean models can lead to a reduction of model biases (e.g., Roberts et al., 2020; Hewitt et al., 2020, 2022; Beech et al., 2022), other issues arise, particularly as to whether mesoscale eddy parameterisations should be dispensed with or employed, and if employed, how they are employed. Without parameterisations, it is known that the eddy processes are mis-represented, and the associated eddy-mean feedback is too weak. On the other hand, it is known that GM-based parameterisations utilised as-is end up damping explicit eddy fluctuations. One suggestion then is to switch off the GM-based parameterisations when the model is regarded as eddy-permitting, via the use of a resolution function (e.g., Hallberg, 2013). Another is to accept there is some damping by GM-based schemes, but re-energise via backscatter approaches (e.g., Jansen et al., 2015b,a, 2019; Bachman, 2019). A proposal considers anisotropic versions of the GM-based parameterisation (e.g., Smith and Gent, 2004). Some works advocate for backscatter only (e.g., Juricke et al., 2020; Chang et al., 2023; Yankovsky et al., 2023), arguing that backscatter re-energises the explicit eddy activity and catalyses for the extra eddy-mean interaction.

In the present work, the principal focus is on quantifying the benefits and/or deficiencies arising from a procedure suggested by Mak et al. (2023) termed field splitting, which modifies the way GM-based schemes are applied. No resolution function, anisotropic GM, backscatter or isoneutral diffusion is employed in this work, although the splitting procedure is not necessarily mutually exclusive of those modelling choices. We apply the splitting procedure to the model and experimental procedure of Ruan et al. (2023) in the eddy-permitting regime, where we expect improvements for reasons to be detailed.

The present work further tests the splitting procedure for a model in a different ocean-relevant setting (a representative mid-latitude gyre system), and comprehensively assesses both the modelled physical and biogeochemical responses, to inform future works using more realistic physical settings, biogeochemical models, and/or forcing scenarios of more direct relevance to climate projection exercises.

In Section 2, we provide the technical problem statement relating to the use of a GM-based parameterisation in an eddy-permitting regime, as well as the motivation and overview for the splitting procedure of Mak et al. (2023). We leverage the model and experimental procedure based largely on the previous work of Ruan et al. (2023) to test our scientific hypotheses, but with the caveat that the results we present in the main body of the work do not employ any isoneutral diffusion, for reasons to be elaborated upon. In Section 3 we provide a detailed comparison of the physical and biogeochemical responses in the set of models under a control pre-industrial scenario. In Section 4 we consider the analogous responses under an idealised climate change scenario, and additionally assess the associated sensitivities. The paper closes with conclusions and discussions in Section 5, detailing implications, limitations and outlooks in light of the present results. In Appendix, we reiterate and further elaborate on why the present work does not include isoneutral diffusion, and provide sample numerical evidence on why there are subtleties with the use of isoneutral diffusion together with state-aware GM-based parameterisations, such as that employed in this work.

## 2. Problem statement, field splitting approach, and numerical set up

### 2.1. Problem statement

The underlying problem in this work relates to whether and how a GM-based parameterisation should be used when the model allows for an explicit representation of mesoscale eddies. Consider a Reynolds decomposition of the velocity (specifically, the advective velocity in the tracer equation only)

$$u = \bar{u} + u' + u^*, \quad \bar{u}' = 0, \qquad (1)$$

where an overbar denotes a Reynolds average (time-average is considered in this work), a prime denotes a deviation from that average such that the average of the deviation is zero (so $u'$ is associated with the explicit eddies), and a star denotes any parameterised component we may wish to add on (so $u^*$ is associated with the parameterised eddies). In the eddy-rich/resolving case without a parameterisation, we would set $u^* = 0$. In a coarse resolution calculation, $u'$ would be effectively zero, and we might mimic the effect of the missing $u'$ by $u^*$, such as via the GM specification (e.g., Gent et al., 1995)

$$u^* = \nabla \times (e_z \times \kappa_{gm}s), \quad s = \frac{\nabla_H \rho}{-\partial \rho / \partial z}. \qquad (2)$$

where $\nabla_H$ is the gradient operator ($\nabla_H$ is the gradient operator only in the horizontal), $e_z$ is the unit vector pointing in the vertical, $\kappa_{gm}$ is the GM or eddy-induced velocity coefficient (in units of $m^2 s^{-1}$), $s$ is the vector encoding the isopycnal slopes in the horizontal directions, and $\rho$ denotes the relevant density variable of interest.

A modelling problem in the eddy-permitting regime is that $u'$, while weak, is not necessarily negligible, and the question is whether $u^*$ should be included or not. If $u^* = 0$ (by setting $\kappa_{gm} = 0$ for example), then the explicit eddies are too weak, leading to a rather weak feedback onto the mean state. As such, we might expect the resulting baroclinic mean flow to be too strong, associated for example with an overly deep stratification (e.g., Fig. 7c of Mak et al., 2023), resulting in an overly strong meridional overturning circulation. A larger lateral transport in the present double-gyre system to be investigated would be expected to lead to meridional heat transport with positive biases, while the
larger vertical transport might be expected to lead to larger Net Primary Production (NPP) via increased nutrient upwelling.

One might suspect some degree of $u^*$ would need to be included. However, that presents another set of problems, in that $u^*$ from the GM-based parameterisations tends to dominate and damps the explicit $u'$ that would be permitted by the model at the relevant spatial resolution. The result is often an eddy-permitting resolution calculation that strongly resembles a coarse resolution model (e.g., Fig 2c of Mak et al., 2023), but at a higher computational cost. The physical rationalisation of this effect is summarised in Fig. 3a of Mak et al. (2023): in the eddy-permitting regime, explicit eddies are still in the geostrophic regime, so velocity fluctuations are associated with isopycnal fluctuations via the thermal wind shear relation. Since GM-based schemes act to mimic baroclinic instability by flattening isopycnals, the isopycnal fluctuations associated with explicit eddies are rapidly damped by the GM-based schemes. As a result, we would expect the performance of the double-gyre model to largely mimic that of the coarse resolution model, with an overturning circulation that is too weak, with negative biases in heat transport, nutrient transport and associated NPP (cf. coarse resolution simulations in Ruan et al., 2023).

To combat the overly diffuse nature of the model in the eddy-permitting regime (without or with GM-based schemes active), approaches based on momentum backscatter have been proposed (e.g., Bachman, 2019; Jansen et al., 2019; Juricke et al., 2020), with the idea that backscatter would strengthen the explicit eddies and catalyse for the extra eddy-mean feedbacks (e.g., Chang et al., 2023; Yankovsky et al., 2023). While there are model improvements as a result of employing only backscatter, it remains to be convincingly demonstrated that backscatter approaches really are supplementing for the extra eddy-mean interaction in the intended fashion.

## 2.2. A field splitting approach for eddy-permitting calculations

The recent work of Mak et al. (2023) revisits the issue of the co-existence of $u'$ and $u^*$, and argues that $u^*$ from a GM-based scheme arises from the collective action of eddies over some sufficiently large-scale region, acting over that same region that should be regarded as eddy-free after the averaging procedure. Since the use of a GM-based scheme as-is in an eddy-permitting model where the state explicitly represents eddies violates the initial working assumption, it is perhaps not surprising that the GM-based schemes have associated modelling deficiencies in the eddy-permitting regime. If one accepts that argument, then a relatively simple fix would be a field splitting approach. Consider a decomposition of the density field as

$$ \rho = \rho_L + \rho_S, \quad (3) $$

where $\rho_L$ is some large-scale density field associated with some non-eddying field, and $\rho_S$ is the residual between the full and large-scale density field. We simply use $\rho_L$ as the input field for the parameterisation (i.e., let $\rho \to \rho_L$ and $s \to s_L$ in Eq. (2), which gives $u_L^* = \nabla \times (e_z \times \kappa_{gm}s_L)$, and proceed as usual. The observation is that $u_L^*$ is now fundamentally a large-scale object at a smaller magnitude, because the input field is large-scale, and the gradient of a quantity that is smooth on the large-scale is of smaller magnitude at the grid-scale. Such an approach should in principle allow $u'$ and $u_L^*$ to co-exist (see the schematic given in Fig. 3b of Mak et al., 2023).

The work of Mak et al. (2023) reports that, in an idealised Southern Ocean configuration, the mean-state improvements appear to require the use of the splitting procedure, and the improvements cannot be attained by switching off the parameterisation or by tuning the parameterisation in the absence of splitting (cf. Fig. 7 Mak et al., 2023). In the double-gyre model at the eddy-permitting resolution to be detailed, relative to the model without a GM-based parameterisation, we would expect that the use of splitting together with a GM-based parameterisation reduces the overly deep stratification, the overly strong meridional overturning circulation, and the overly large NPP resulting

from the larger nutrient transport. Further, we hypothesise that the use of splitting would lead to improvements that cannot be achieved by a model with a GM-based parameterisation without splitting via tuning the free parameters.

## 2.3. Model set up

To test for the aforementioned hypotheses and quantify the benefits and/or deficiencies afforded to the modelled physical and biogeochemical responses by the splitting procedure, we employ a double-gyre model coupled to an idealised biogeochemistry model LOBSTER (cf. Lévy et al., 2010, 2012); the model is essentially that reported in Couespel et al. (2021) and Ruan et al. (2023), and here we only recap the essentials (see Ruan et al., 2023, for in-depth details). The model is created in the Nucleus for European Modelling of the Ocean framework (NEMO, version 4.0.5 r14538; Madec 2008). The domain of the double-gyre model is a ‘straightened’ version of the NEMO rotated gyre configuration test case. The square domain has sides of length $L = 3180$ km, employing a regular horizontal grid-spacing. The model is on a $\beta$-plane centred at around 35° N, extending to 20° N to the south and 50° N to the north. The lateral boundary momentum condition of each side is set to be free slip, and a non-linear drag is applied on the bottom boundary.

Given that the bathymetry can have impacts on the large-scale circulation and thus affect the physical and biogeochemical representations, we consider a model with a slope on the west and east sides in this study, in contrast to the previous works of Couespel et al. (2021) and Ruan et al. (2023) with no bathymetry. Following the work of Jackson et al. (2006), we consider a bottom bathymetry that varies only in the zonal direction, given by (in units of m)

$$ H = 100 + (H_i - 100) \left[ 1 - e^{-x^2/\sigma^2} - e^{-(L-x)^2/\sigma^2} \right], \quad (4) $$

where $H_i = 4000$ m is the total depth of the modelled ocean (the ocean is shallowest at 100 m), $x$ is the offshore distance from the western boundary, $L$ is the width of the domain, and $\sigma$ is taken to be 100,000 m = 100 km (roughly corresponding to the choice of 1° employed in Jackson et al. 2006). In the NEMO model, 31 uneven vertical layers are employed, in line with the standard gyre test case, and the vertical coordinate is chosen to be in z-coordinate with partial steps. The main consequence of employing such a bathymetry is a slight southward shift of the modelled Western Boundary Current relative to the relevant previous works with no bathymetry (Couespel et al., 2021; Ruan et al., 2023). The qualitative conclusions we draw from this work are robust and carry over to the case with no bathymetry (not shown).

The key model parameter settings of the set of calculations are displayed in Table 1. Tracer and momentum lateral diffusion are in the geopotential direction for all models reported here. Tracer advection is still with the MUSCL scheme. Momentum advection is processed in vector form with a second order centred scheme, and the vertical mixing is parameterised by the turbulent kinetic energy closure of Gaspar et al. (1990). We keep the choice of linear equation of state, so that the splitting approach as implemented in Mak et al. (2023) using a filtered temperature and salinity field to compute the resulting large-scale density field $\rho_L$ can be used without further modifications. The zonally symmetric atmospheric forcing employs the flux formulation, with a predetermined repeating seasonal cycle for the wind stress, penetrative solar radiation, pseudo-atmospheric temperature and freshwater flux. For more details, see Fig. 1 of Lévy et al. (2010) or Ruan et al. (2023).

To reiterate, no isoneutral diffusion is employed in the results provided in the main body of the present article. While this might be counter to standard practices, the reason is that there appears to be a positive feedback loop when isoneutral diffusion is utilised with state-aware GM-based parameterisations, leading to further changes in the stratification, affecting the various physical and biogeochemical responses. The positive feedback loop is elaborated upon in Appendix,
larger vertical transport might be expected to lead to larger Net Primary Production (NPP) via increased nutrient upwelling.

One might suspect some degree of $u^*$ would need to be included. However, that presents another set of problems, in that $u^*$ from the GM-based parameterisations tends to dominate and damps the explicit $u'$ that would be permitted by the model at the relevant spatial resolution. The result is often an eddy-permitting resolution calculation that strongly resembles a coarse resolution model (e.g., Fig 2c of Mak et al., 2023), but at a higher computational cost. The physical rationalisation of this effect is summarised in Fig. 3a of Mak et al. (2023): in the eddy-permitting regime, explicit eddies are still in the geostrophic regime, so velocity fluctuations are associated with isopycnal fluctuations via the thermal wind shear relation. Since GM-based schemes act to mimic baroclinic instability by flattening isopycnals, the isopycnal fluctuations associated with explicit eddies are rapidly damped by the GM-based schemes. As a result, we would expect the performance of the double-gyre model to largely mimic that of the coarse resolution model, with an overturning circulation that is too weak, with negative biases in heat transport, nutrient transport and associated NPP (cf. coarse resolution simulations in Ruan et al., 2023).

To combat the overly diffuse nature of the model in the eddy-permitting regime (without or with GM-based schemes active), approaches based on momentum backscatter have been proposed (e.g., Bachman, 2019; Jansen et al., 2019; Juricke et al., 2020), with the idea that backscatter would strengthen the explicit eddies and catalyse for the extra eddy-mean feedbacks (e.g., Chang et al., 2023; Yankovsky et al., 2023). While there are model improvements as a result of employing only backscatter, it remains to be convincingly demonstrated that backscatter approaches really are supplementing for the extra eddy-mean interaction in the intended fashion.

## 2.2. A field splitting approach for eddy-permitting calculations

The recent work of Mak et al. (2023) revisits the issue of the co-existence of $u'$ and $u^*$, and argues that $u^*$ from a GM-based scheme arises from the collective action of eddies over some sufficiently large-scale region, acting over that same region that should be regarded as eddy-free after the averaging procedure. Since the use of a GM-based scheme as-is in an eddy-permitting model where the state explicitly represents eddies violates the initial working assumption, it is perhaps not surprising that the GM-based schemes have associated modelling deficiencies in the eddy-permitting regime. If one accepts that argument, then a relatively simple fix would be a field splitting approach. Consider a decomposition of the density field as

$$ \rho = \rho_L + \rho_S, \qquad (3) $$

where $\rho_L$ is some large-scale density field associated with some non-eddying field, and $\rho_S$ is the residual between the full and large-scale density field. We simply use $\rho_L$ as the input field for the parameterisation (i.e., let $\rho \to \rho_L$ and $s \to s_L$ in Eq. (2), which gives $u_L^* = \nabla \times (e_z \times \kappa_{gm}s_L)$, and proceed as usual. The observation is that $u_L^*$ is now fundamentally a large-scale object at a smaller magnitude, because the input field is large-scale, and the gradient of a quantity that is smooth on the large-scale is of smaller magnitude at the grid-scale. Such an approach should in principle allow $u'$ and $u_L^*$ to co-exist (see the schematic given in Fig. 3b of Mak et al., 2023).

The work of Mak et al. (2023) reports that, in an idealised Southern Ocean configuration, the mean-state improvements appear to require the use of the splitting procedure, and the improvements cannot be attained by switching off the parameterisation or by tuning the parameterisation in the absence of splitting (cf. Fig. 7 Mak et al., 2023). In the double-gyre model at the eddy-permitting resolution to be detailed, relative to the model without a GM-based parameterisation, we would expect that the use of splitting together with a GM-based parameterisation reduces the overly deep stratification, the overly strong meridional overturning circulation, and the overly large NPP resulting

from the larger nutrient transport. Further, we hypothesise that the use of splitting would lead to improvements that cannot be achieved by a model with a GM-based parameterisation without splitting via tuning the free parameters.

## 2.3. Model set up

To test for the aforementioned hypotheses and quantify the benefits and/or deficiencies afforded to the modelled physical and biogeochemical responses by the splitting procedure, we employ a double-gyre model coupled to an idealised biogeochemistry model LOBSTER (cf. Lévy et al., 2010, 2012); the model is essentially that reported in Couespel et al. (2021) and Ruan et al. (2023), and here we only recap the essentials (see Ruan et al., 2023, for in-depth details). The model is created in the Nucleus for European Modelling of the Ocean framework (NEMO, version 4.0.5 r14538; Madec 2008). The domain of the double-gyre model is a ‘straightened’ version of the NEMO rotated gyre configuration test case. The square domain has sides of length $L = 3180$ km, employing a regular horizontal grid-spacing. The model is on a $\beta$-plane centred at around 35° N, extending to 20° N to the south and 50° N to the north. The lateral boundary momentum condition of each side is set to be free slip, and a non-linear drag is applied on the bottom boundary.

Given that the bathymetry can have impacts on the large-scale circulation and thus affect the physical and biogeochemical representations, we consider a model with a slope on the west and east sides in this study, in contrast to the previous works of Couespel et al. (2021) and Ruan et al. (2023) with no bathymetry. Following the work of Jackson et al. (2006), we consider a bottom bathymetry that varies only in the zonal direction, given by (in units of m)

$$ H = 100 + (H_i - 100) \left[ 1 - e^{-x^2/\sigma^2} - e^{-(L-x)^2/\sigma^2} \right], \qquad (4) $$

where $H_i = 4000$ m is the total depth of the modelled ocean (the ocean is shallowest at 100 m), $x$ is the offshore distance from the western boundary, $L$ is the width of the domain, and $\sigma$ is taken to be 100,000 m = 100 km (roughly corresponding to the choice of 1° employed in Jackson et al. 2006). In the NEMO model, 31 uneven vertical layers are employed, in line with the standard gyre test case, and the vertical coordinate is chosen to be in z-coordinate with partial steps. The main consequence of employing such a bathymetry is a slight southward shift of the modelled Western Boundary Current relative to the relevant previous works with no bathymetry (Couespel et al., 2021; Ruan et al., 2023). The qualitative conclusions we draw from this work are robust and carry over to the case with no bathymetry (not shown).

The key model parameter settings of the set of calculations are displayed in Table 1. Tracer and momentum lateral diffusion are in the geopotential direction for all models reported here. Tracer advection is still with the MUSCL scheme. Momentum advection is processed in vector form with a second order centred scheme, and the vertical mixing is parameterised by the turbulent kinetic energy closure of Gaspar et al. (1990). We keep the choice of linear equation of state, so that the splitting approach as implemented in Mak et al. (2023) using a filtered temperature and salinity field to compute the resulting large-scale density field $\rho_L$ can be used without further modifications. The zonally symmetric atmospheric forcing employs the flux formulation, with a predetermined repeating seasonal cycle for the wind stress, penetrative solar radiation, pseudo-atmospheric temperature and freshwater flux. For more details, see Fig. 1 of Lévy et al. (2010) or Ruan et al. (2023).

To reiterate, no isoneutral diffusion is employed in the results provided in the main body of the present article. While this might be counter to standard practices, the reason is that there appears to be a positive feedback loop when isoneutral diffusion is utilised with state-aware GM-based parameterisations, leading to further changes in the stratification, affecting the various physical and biogeochemical responses. The positive feedback loop is elaborated upon in Appendix,
<table><caption>**Table 1**
Key model parameter values of the set of calculations.</caption><thead><tr><th></th><th>1/4° (R4, SPLIT and GEOM)</th><th>1/12° (R12)</th></tr></thead><tbody><tr><td>Horizontal resolution Δx = Δy</td><td>26.5 km</td><td>8.83 km</td></tr><tr><td>Time step Δt</td><td>20 mins</td><td>10 mins</td></tr><tr><td>Tracer diffusion κ<sub>T</sub></td><td>Geopotential bi-Laplacian ∇<sup>4</sup>, -5 × 10<sup>9</sup> m s<sup>-4</sup></td><td>Geopotential bi-Laplacian ∇<sup>4</sup>, -10<sup>9</sup> m s<sup>-4</sup></td></tr><tr><td>Momentum diffusion κ<sub>M</sub></td><td>Geopotential bi-Laplacian ∇<sup>4</sup>, -2.5 × 10<sup>11</sup> m s<sup>-4</sup></td><td>Geopotential bi-Laplacian ∇<sup>4</sup>, -3 × 10<sup>10</sup> m s<sup>-4</sup></td></tr><tr><td>Tracer advection</td><td>MUSCL</td><td>MUSCL</td></tr><tr><td>Momentum advection</td><td>Centerd standard kinetic energy scheme</td><td>Centerd standard kinetic energy scheme</td></tr></tbody></table>

<table><thead><tr><th></th><th>SPLIT</th><th>GEOM</th></tr></thead><tbody><tr><td>Tuning parameter α</td><td>0.065</td><td>0.065</td></tr><tr><td>Energy diffusivity η<sub>E</sub></td><td>500 m<sup>2</sup> s<sup>-1</sup></td><td>500 m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Minimum energy level E<sub>0</sub></td><td>1.0 m<sup>3</sup> s<sup>-2</sup></td><td>1.0 m<sup>3</sup> s<sup>-2</sup></td></tr><tr><td>Dissipation time-scale λ<sup>-1</sup></td><td>135 days</td><td>135 days</td></tr><tr><td>Filter length scale L</td><td>100 km</td><td>—</td></tr><tr><td>Pre-conditioning param. γ</td><td>75</td><td>—</td></tr></tbody></table>

where we also provide supporting numerical evidence for such a positive feedback loop, but also isoneutral diffusion can lead to non-negligible damping of the explicit variability in the eddy-permitting regime. Both effects lead to non-trivial model responses, and we focus only on the GM-based parameterisation in the present article for simplicity.

The model employs the simplified biogeochemistry model LOBSTER, which takes nitrogen as the currency, solving six biogeochemical variables of phytoplankton, zooplankton, detritus, dissolved organic matter, nitrate and ammonium (e.g., Lévy et al., 2012). The uptake of nitrate and ammonium by phytoplankton determined as the Net Primary Production (NPP) is of particular interest in this study. The absence of physiological changes represented in the idealised model means the plankton is primarily affected by the modelled flow and the related advective tendencies, and a dominant control on NPP is the nutrient supply over the large length-scales via the modelled circulation.

For the present work, the GM-based parameterisation that we primarily focus on is the GM-version of the GEOMETRIC (Marshall et al., 2012; Mak et al., 2018, 2022b), given its use in the previous work of Ruan et al. (2023) for the coarse resolution case, and other works that have demonstrated the use of GEOMETRIC leads to various improvements in the modelled mean state (Mak et al., 2018, 2022b, 2023; Wei et al., 2024). Via the analysis given in Mak et al. (2023), if the splitting approach (given essentially by Eq. (2) with $ρ → ρ_L$ and $s → s_L$) is to be used with GEOMETRIC, for consistency we should take

$$ \kappa_{gm} = \alpha \frac{\int E \, dz}{\int (M_L^2/N_L) \, dz}, \qquad (5) $$

where $α$ is a non-dimensional tuning parameter ($|α| ≤ 1$), $E$ is the total (potential and kinetic) parameterised eddy energy, $M_L^2 ∼ |∇_H ρ_L|$ is the large-scale horizontal buoyancy frequency, and $N_L^2 ∼ -∂ρ_L/∂z$ is the large-scale vertical buoyancy frequency. The vertical integration in Eq. (5) results in a time and horizontally varying $κ_{gm}(x, y, t)$, and the depth-integrated total parameterised eddy energy is prognostically constrained by the eddy energy budget

$$ \frac{d}{dt} \int E \, dz + \nabla_H \cdot \left( (\tilde{u}^z - |c| \mathbf{e}_x) \int E \, dz \right) = \underbrace{\int \kappa_{gm} \frac{M_L^4}{N_L^2} dz - \lambda \int (E - E_0) dz + \eta_E \nabla_H^2 \int E \, dz}_{\text{source}} \qquad (6) $$

Here, the source of eddy energy with the use of a splitting approach stems from large-scale slumping of density surfaces. The depth-integrated eddy energy is still advected by the depth average flow

$\tilde{u}^z$, with westward propagation at the long Rossby wave phase speed $|c|$ (e.g., Chelton et al., 2011; Klockner and Marshall, 2014), diffused in the horizontal (Grooms, 2015; Ni et al., 2020a,b) with eddy energy diffusivity $η_E$, dissipated at the rate $λ$ ($λ^{-1}$ is an eddy energy dissipation time-scale, which can in principle vary in time and space; cf. Mak et al. 2022a; Torres et al. 2023; Wilder et al. 2023), and $E_0$ is a minimum eddy energy level. The parameterisation appears to be scale-aware in the eddy energetics when splitting is employed, since the co-existing explicit eddy feedback affects the resolved mean-state, and the resolved mean-state modifies the parameterised eddy feedback via changes in the parameterised eddy energetics (Mak et al., 2023).

## 2.4. Experimental set up

The model spin up follows the same procedure as that detailed in Ruan et al. (2023). The initial spun-up state from a 1° resolution with a constant $κ_{gm} = 1000 \, m^2 \, s^{-1}$ and a flat bottom, starting from model year -2300 to year -300 that already exists from Ruan et al. (2023) is interpolated onto the new domain and appropriately masked, and at model year -300 four sets of perturbation experiments were considered, running up to model year 0. We have considered sample calculations where we spin-up from rest on the sloped domain; while there are differences in the deeper parts of the ocean, over the top 700 m where we compute our bulk diagnostics the differences are minor (not shown). For the perturbation experiments, a 1/12° horizontal resolution model (R12) resolving most of the mesoscale eddies serves as a model truth as a reference. The primary focus here are the three eddy-permitting calculations at 1/4° horizontal resolution:

* R4, with no GM-based parameterisation active,

* SPLIT, with the GM-based GEOMETRIC parameterisation and with splitting (i.e., Eqs. (5) and (6)),

* GEOM, with the GM-based GEOMETRIC parameterisation but with no splitting (i.e., Eqs. (5) and (6) without the subscript $L$; cf. Eq. 3–4 in Ruan et al. 2023).

For each of the experiments, a pre-industrial control scenario (assigned a suffix CTL) and an idealised climate change scenario (assigned a suffix CC) are performed from model year 0 to 70. The pre-industrial control takes the standard forcing as-is, while in the idealised climate change scenario the atmospheric pseudo-temperature has an added linear trend of +0.04 °C yr⁻¹ over the 70 model years, to mimic the SSP5-8.5 scenario in the North Atlantic (e.g., Tokarska et al., 2020). All time-averaged diagnostics reported in this work are based on the last five years of the simulation (spanning from the start of model year 66 and the end of model year 70).

The GEOMETRIC parameters are given in Table 1 and are the same in the SPLIT and GEOM calculations, with $α$ chosen to be closer
<table><caption>**Table 1**
Key model parameter values of the set of calculations.</caption><thead><tr><th></th><th>1/4° (R4, SPLIT and GEOM)</th><th>1/12° (R12)</th></tr></thead><tbody><tr><td>Horizontal resolution Δx = Δy</td><td>26.5 km</td><td>8.83 km</td></tr><tr><td>Time step Δt</td><td>20 mins</td><td>10 mins</td></tr><tr><td>Tracer diffusion κ<sub>T</sub></td><td>Geopotential bi-Laplacian ∇<sup>4</sup>, -5 × 10<sup>9</sup> m s<sup>-4</sup></td><td>Geopotential bi-Laplacian ∇<sup>4</sup>, -10<sup>9</sup> m s<sup>-4</sup></td></tr><tr><td>Momentum diffusion κ<sub>M</sub></td><td>Geopotential bi-Laplacian ∇<sup>4</sup>, -2.5 × 10<sup>11</sup> m s<sup>-4</sup></td><td>Geopotential bi-Laplacian ∇<sup>4</sup>, -3 × 10<sup>10</sup> m s<sup>-4</sup></td></tr><tr><td>Tracer advection</td><td>MUSCL</td><td>MUSCL</td></tr><tr><td>Momentum advection</td><td>Centerd standard kinetic energy scheme</td><td>Centerd standard kinetic energy scheme</td></tr></tbody></table>

<table><thead><tr><th></th><th>SPLIT</th><th>GEOM</th></tr></thead><tbody><tr><td>Tuning parameter α</td><td>0.065</td><td>0.065</td></tr><tr><td>Energy diffusivity η<sub>E</sub></td><td>500 m<sup>2</sup> s<sup>-1</sup></td><td>500 m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Minimum energy level E<sub>0</sub></td><td>1.0 m<sup>3</sup> s<sup>-2</sup></td><td>1.0 m<sup>3</sup> s<sup>-2</sup></td></tr><tr><td>Dissipation time-scale λ<sup>-1</sup></td><td>135 days</td><td>135 days</td></tr><tr><td>Filter length scale L</td><td>100 km</td><td>–</td></tr><tr><td>Pre-conditioning param. γ</td><td>75</td><td>–</td></tr></tbody></table>

where we also provide supporting numerical evidence for such a positive feedback loop, but also isoneutral diffusion can lead to non-negligible damping of the explicit variability in the eddy-permitting regime. Both effects lead to non-trivial model responses, and we focus only on the GM-based parameterisation in the present article for simplicity.

The model employs the simplified biogeochemistry model LOBSTER, which takes nitrogen as the currency, solving six biogeochemical variables of phytoplankton, zooplankton, detritus, dissolved organic matter, nitrate and ammonium (e.g., Lévy et al., 2012). The uptake of nitrate and ammonium by phytoplankton determined as the Net Primary Production (NPP) is of particular interest in this study. The absence of physiological changes represented in the idealised model means the plankton is primarily affected by the modelled flow and the related advective tendencies, and a dominant control on NPP is the nutrient supply over the large length-scales via the modelled circulation.

For the present work, the GM-based parameterisation that we primarily focus on is the GM-version of the GEOMETRIC (Marshall et al., 2012; Mak et al., 2018, 2022b), given its use in the previous work of Ruan et al. (2023) for the coarse resolution case, and other works that have demonstrated the use of GEOMETRIC leads to various improvements in the modelled mean state (Mak et al., 2018, 2022b, 2023; Wei et al., 2024). Via the analysis given in Mak et al. (2023), if the splitting approach (given essentially by Eq. (2) with $ρ → ρ_L$ and $s → s_L$) is to be used with GEOMETRIC, for consistency we should take

$$ \kappa_{gm} = \alpha \frac{\int E \, dz}{\int (M_L^2/N_L) \, dz}, \qquad (5) $$

where $α$ is a non-dimensional tuning parameter ($|α| ≤ 1$), $E$ is the total (potential and kinetic) parameterised eddy energy, $M_L^2 ∼ |∇_H ρ_L|$ is the large-scale horizontal buoyancy frequency, and $N_L^2 ∼ -∂ρ_L/∂z$ is the large-scale vertical buoyancy frequency. The vertical integration in Eq. (5) results in a time and horizontally varying $κ_{gm}(x, y, t)$, and the depth-integrated total parameterised eddy energy is prognostically constrained by the eddy energy budget

$$ \frac{d}{dt} \int E \, dz + \nabla_H \cdot \left( (\tilde{u}^z - |c| \mathbf{e}_x) \int E \, dz \right) = \underbrace{\int \kappa_{gm} \frac{M_L^4}{N_L^2} dz - \lambda \int (E - E_0) dz + \eta_E \nabla_H^2 \int E \, dz}_{\text{source}} \qquad (6) $$

Here, the source of eddy energy with the use of a splitting approach stems from large-scale slumping of density surfaces. The depth-integrated eddy energy is still advected by the depth average flow

$\tilde{u}^z$, with westward propagation at the long Rossby wave phase speed $|c|$ (e.g., Chelton et al., 2011; Klockner and Marshall, 2014), diffused in the horizontal (Grooms, 2015; Ni et al., 2020a,b) with eddy energy diffusivity $η_E$, dissipated at the rate $λ$ ($λ^{-1}$ is an eddy energy dissipation time-scale, which can in principle vary in time and space; cf. Mak et al. 2022a; Torres et al. 2023; Wilder et al. 2023), and $E_0$ is a minimum eddy energy level. The parameterisation appears to be scale-aware in the eddy energetics when splitting is employed, since the co-existing explicit eddy feedback affects the resolved mean-state, and the resolved mean-state modifies the parameterised eddy feedback via changes in the parameterised eddy energetics (Mak et al., 2023).

## 2.4. Experimental set up

The model spin up follows the same procedure as that detailed in Ruan et al. (2023). The initial spun-up state from a 1° resolution with a constant $κ_{gm} = 1000 \, m^2 \, s^{-1}$ and a flat bottom, starting from model year -2300 to year -300 that already exists from Ruan et al. (2023) is interpolated onto the new domain and appropriately masked, and at model year -300 four sets of perturbation experiments were considered, running up to model year 0. We have considered sample calculations where we spin-up from rest on the sloped domain; while there are differences in the deeper parts of the ocean, over the top 700 m where we compute our bulk diagnostics the differences are minor (not shown). For the perturbation experiments, a 1/12° horizontal resolution model (R12) resolving most of the mesoscale eddies serves as a model truth as a reference. The primary focus here are the three eddy-permitting calculations at 1/4° horizontal resolution:

* R4, with no GM-based parameterisation active,

* SPLIT, with the GM-based GEOMETRIC parameterisation and with splitting (i.e., Eqs. (5) and (6)),

* GEOM, with the GM-based GEOMETRIC parameterisation but with no splitting (i.e., Eqs. (5) and (6) without the subscript $L$; cf. Eq. 3–4 in Ruan et al. 2023).

For each of the experiments, a pre-industrial control scenario (assigned a suffix CTL) and an idealised climate change scenario (assigned a suffix CC) are performed from model year 0 to 70. The pre-industrial control takes the standard forcing as-is, while in the idealised climate change scenario the atmospheric pseudo-temperature has an added linear trend of +0.04 °C yr⁻¹ over the 70 model years, to mimic the SSP5-8.5 scenario in the North Atlantic (e.g., Tokarska et al., 2020). All time-averaged diagnostics reported in this work are based on the last five years of the simulation (spanning from the start of model year 66 and the end of model year 70).

The GEOMETRIC parameters are given in Table 1 and are the same in the SPLIT and GEOM calculations, with $α$ chosen to be closer
<figure><img src="image_2.png" /><figcaption>Fig. 1. (a) Model bathymetry (shading) and the large-scale potential vorticity contours as f/H (m⁻¹s⁻¹) in the present model, where f is the Coriolis parameter and H is the depth. The diagnosed surface temperature field (°C) from SPLIT, showing (b) the full temperature field, (c) the associated large-scale temperature field, and (d) residual or small-scale temperature field.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. The resulting κ gm (x, y) (m²s⁻¹) distribution from (a) SPLIT and (b) GEOM, with some representative contours of κ gm marked on.</figcaption></figure>

to the recent works employing GEOMETRIC (Ruan et al., 2023; Mak et al., 2023; Wei et al., 2024). The splitting procedure, where active, is performed every model day using the procedure detailed in Mak et al. (2023): briefly, a diffusion based filter in space with a pre-defined length-scale L (taken to be 100 km here) is performed per model level to obtain ρ$_L$ , where the filtering kernel is closely related to the Matérn auto-covariance (e.g., Whittle, 1963; Lindgren et al., 2018). The resulting diffusion problem is solved with a Richardson pre-conditioning with regularisation parameter γ = 75 every model day. Fig. 1a shows the bathymetry and the resulting f/H contours (that the geostrophic flow should be constrained somewhat to follow), noting the poleward deflection of the contours as we move eastwards from the western boundary. A sample output from SPLIT is given in Fig. 1b- d for a snapshot of the sea surface temperature, showing the total, filtered and residual field respectively, demonstrating that the splitting procedure leaves a portion of the explicit fluctuations intact, as seen in the residual field. The analogous results in R12, GEOM and R4 have been omitted here for brevity: R12 and R4 permits explicit fluctuations, while GEOM largely resembles the coarse resolution calculations (cf. Ruan et al., 2023; Mak et al., 2023).

In Fig. 2 we show the resulting κ gm (x, y) field for both SPLIT and GEOM. Both are large in the Western Boundary Current region, because of large simulated total parameterised eddy energy E . Notice that values of κ gm are much more modest in SPLIT than in GEOM (domain- averaged value at 468 and 2361 m² s⁻¹ respectively), given the same parameter choices. While the reduction in the horizontal gradients of the associated filtered density field ρ$_L$ does lead to a reduction of the simulated E (by affecting κ gm ), a fundamental difference leading to non-damping of explicit eddies is that the eddy-induced velocity u ∗ L is a large-scale rather than grid-scale object. The reported behaviour later is found to crucially depend on the use of u ∗ L (i.e., splitting), and less on a reduction in κ gm (which can be achieved by tuning α and/or λ ).

One point we make is that, since we have a non-trivial bathymetry in the present model set up, it might be possible that a tapering of the GM coefficient κ gm is required over continental slopes, where the

dynamics of eddies and topographic effects can differ from the open ocean (e.g., LaCasc e and Brink, 2000; Stewart and Thompson, 2013; Wang and Stewart, 2018). The present implementation of GEOMETRIC in NEMO takes a simple choice of tapering κ gm to zero as the Rossby deformation becomes sufficiently small, which in the present model is largely dictated by the modelled water depth. Simulations with and without tapering in the present model (enabled by commenting out the relevant lines of the source code) seem to make no qualitative difference to any of our reported results. There are more advanced choices based on slope parameters or Burger numbers that have shown promise (e.g., Wei et al., 2022, 2024; Nummelin and Isachsen, 2024), which may affect the model response, although we have not considered implementations of those procedures here. The present reported results have the Rossby number based tapering of κ gm deactivated.

3. Pre-industrial control responses

3.1. Physical responses

We first show some metrics relating to the modelled circulation associated with the set of calculations. In Fig. 3 we show the barotropic streamfunction Ψ$_{baro}$, calculated as

$$
\Psi_{\text{baro}} = \int_{\bar{x}=0}^{\bar{x}=x} \int_{-H}^{0} \bar{v}(\bar{x}, y, z) \, dz \, d\bar{x}. \qquad (7)
$$

The streamfunction Ψ$_{baro}$ displays the basic features of a subtropical gyre to the south and a subpolar gyre to the north, separated by a Western Boundary Current region. Compared to the analogous diag- nostic reported in Fig. 2 in Ruan et al. (2023), a main difference here is in the poleward deflection of the subpolar gyre as we move away from the western boundary, consistent with the presence of the sloping bathymetry (e.g., Jackson et al., 2006). The R12, R4 and SPLIT calculations all show some explicit representation of a Western Boundary Current as well as some semblance of fluctuations even in the time-averaged data, unlike the GEOM calculation, which largely resembles the coarse resolution calculations (cf. Ruan et al., 2023, Fig. 2b). Examination of the eddy kinetic energy field or snapshots of surface relative vorticity field indicates that the GEOM calculations possess very weak fluctuations for the present choice of parameter values, relative to the R4 and SPLIT calculations (not shown). We note that the R12 Western Boundary Current is still slightly south of the latitudinal centre line even though the zonal wind stress is symmetric about the latitude centre line, and extends more eastward (cf. Fig. 3a here, and Fig. 2a of Ruan et al. 2023 for R12). Moreover, all calculations display some representation of the re-circulating Fofonoff gyres towards the northern and southern boundary, although the northern one is somewhat weaker, presumably due to the presence of the non-trivial bathymetry in this work.

A subtlety in this work is that the eddy-permitting calculations do have a representation of the Western Boundary Current in some way (at least for R4 and SPLIT), affecting the size of what would be identified as
<figure><img src="image_3.png" /><figcaption>Fig. 1. (a) Model bathymetry (shading) and the large-scale potential vorticity contours as f/H (m⁻¹s⁻¹) in the present model, where f is the Coriolis parameter and H is the depth. The diagnosed surface temperature field (°C) from SPLIT, showing (b) the full temperature field, (c) the associated large-scale temperature field, and (d) residual or small-scale temperature field.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2. The resulting κ gm (x, y) (m²s⁻¹) distribution from (a) SPLIT and (b) GEOM, with some representative contours of κ gm marked on.</figcaption></figure>

to the recent works employing GEOMETRIC (Ruan et al., 2023; Mak et al., 2023; Wei et al., 2024). The splitting procedure, where active, is performed every model day using the procedure detailed in Mak et al. (2023): briefly, a diffusion based filter in space with a pre-defined length-scale L (taken to be 100 km here) is performed per model level to obtain ρ$_L$, where the filtering kernel is closely related to the Matérn auto-covariance (e.g., Whittle, 1963; Lindgren et al., 2018). The resulting diffusion problem is solved with a Richardson pre-conditioning with regularisation parameter γ = 75 every model day. Fig. 1a shows the bathymetry and the resulting f/H contours (that the geostrophic flow should be constrained somewhat to follow), noting the poleward deflection of the contours as we move eastwards from the western boundary. A sample output from SPLIT is given in Fig. 1b-d for a snapshot of the sea surface temperature, showing the total, filtered and residual field respectively, demonstrating that the splitting procedure leaves a portion of the explicit fluctuations intact, as seen in the residual field. The analogous results in R12, GEOM and R4 have been omitted here for brevity: R12 and R4 permits explicit fluctuations, while GEOM largely resembles the coarse resolution calculations (cf. Ruan et al., 2023; Mak et al., 2023).

In Fig. 2 we show the resulting κ gm (x, y) field for both SPLIT and GEOM. Both are large in the Western Boundary Current region, because of large simulated total parameterised eddy energy E. Notice that values of κ gm are much more modest in SPLIT than in GEOM (domain-averaged value at 468 and 2361 m² s⁻¹ respectively), given the same parameter choices. While the reduction in the horizontal gradients of the associated filtered density field ρ$_L$ does lead to a reduction of the simulated E (by affecting κ gm), a fundamental difference leading to non-damping of explicit eddies is that the eddy-induced velocity u* is a large-scale rather than grid-scale object. The reported behaviour later is found to crucially depend on the use of u* (i.e., splitting), and less on a reduction in κ gm (which can be achieved by tuning α and/or λ).

One point we make is that, since we have a non-trivial bathymetry in the present model set up, it might be possible that a tapering of the GM coefficient κ gm is required over continental slopes, where the

dynamics of eddies and topographic effects can differ from the open ocean (e.g., LaCascione and Brink, 2000; Stewart and Thompson, 2013; Wang and Stewart, 2018). The present implementation of GEOMETRIC in NEMO takes a simple choice of tapering κ gm to zero as the Rossby deformation becomes sufficiently small, which in the present model is largely dictated by the modelled water depth. Simulations with and without tapering in the present model (enabled by commenting out the relevant lines of the source code) seem to make no qualitative difference to any of our reported results. There are more advanced choices based on slope parameters or Burger numbers that have shown promise (e.g., Wei et al., 2022, 2024; Nummelin and Isachsen, 2024), which may affect the model response, although we have not considered implementations of those procedures here. The present reported results have the Rossby number based tapering of κ gm deactivated.

## 3. Pre-industrial control responses

### 3.1. Physical responses

We first show some metrics relating to the modelled circulation associated with the set of calculations. In Fig. 3 we show the barotropic streamfunction Ψ$_{baro}$, calculated as

$$ \Psi_{\text{baro}} = \int_{\bar{x}=0}^{\bar{x}=x} \int_{-H}^{0} \bar{v}(\bar{x}, y, z) \, dz \, d\bar{x}. \qquad (7) $$

The streamfunction Ψ$_{baro}$ displays the basic features of a subtropical gyre to the south and a subpolar gyre to the north, separated by a Western Boundary Current region. Compared to the analogous diagnostic reported in Fig. 2 in Ruan et al. (2023), a main difference here is in the poleward deflection of the subpolar gyre as we move away from the western boundary, consistent with the presence of the sloping bathymetry (e.g., Jackson et al., 2006). The R12, R4 and SPLIT calculations all show some explicit representation of a Western Boundary Current as well as some semblance of fluctuations even in the time-averaged data, unlike the GEOM calculation, which largely resembles the coarse resolution calculations (cf. Ruan et al., 2023, Fig. 2b). Examination of the eddy kinetic energy field or snapshots of surface relative vorticity field indicates that the GEOM calculations possess very weak fluctuations for the present choice of parameter values, relative to the R4 and SPLIT calculations (not shown). We note that the R12 Western Boundary Current is still slightly south of the latitudinal centre line even though the zonal wind stress is symmetric about the latitude centre line, and extends more eastward (cf. Fig. 3a here, and Fig. 2a of Ruan et al. 2023 for R12). Moreover, all calculations display some representation of the re-circulating Fofonoff gyres towards the northern and southern boundary, although the northern one is somewhat weaker, presumably due to the presence of the non-trivial bathymetry in this work.

A subtlety in this work is that the eddy-permitting calculations do have a representation of the Western Boundary Current in some way (at least for R4 and SPLIT), affecting the size of what would be identified as
<figure><img src="image_4.png" /><figcaption>Fig. 3. The barotropic streamfunction $\Psi_{\text{baro}}$ (in Sv = $10^6 \text{ m}^3 \text{ s}^{-1}$) of (a) R12, (b) R4, (c) SPLIT, and (d) GEOM, with the zero contour added. The black dashed lines represent the sample latitudinal cross sections $S_{1,2,3}$ going from south to north.</figcaption></figure>

the subpolar gyre. While the works of Couespel et al. (2021) and Ruan et al. (2023) focus on the analysis within the subpolar gyre (defined over fixed geographical locations) primarily because it is the most bio-active region in the modelled domain, in this work we consider the whole domain but exclude the Fofonoff gyres as the analysis region when computing averaged/integrated quantities. For computing fluxes, we consider the cross sections at $y = 25^\circ$ and $45^\circ$ N marked on as black dashed lines in the figures as appropriate (denoted $S_{1,3}$), roughly as the southern boundary of the subtropical gyre, and the northern boundary of the subpolar gyre respectively. For completeness, we also mark on and compute fluxes over the $y = 38^\circ$ N section (denoted $S_2$), which is an empirically determined location that is sufficiently north of the simulated Western Boundary Current (see Fig. 4a-d). The analysis region of primary interest is the region bounded in the horizontal by $S_1$ and $S_3$, and in the vertical by $z = -700 \text{ m}$; we neglect the deeper parts of the ocean since these regions would presumably not have equilibrated within the 70 model year period we are considering.

As another measure of the circulation, we show in Fig. 4 the residual meridional overturning circulation (MOC) streamfunction $\Psi_{\text{MOC}}$, diagnosed as

$$ \Psi_{\text{MOC}} = \int_{z=-H}^{z=H} \int_{0}^{L_x} [\bar{v}(x, y, \bar{z}) + v^*(x, y, \bar{z})] \, dx \, d\bar{z}, \quad (8) $$

as well as a histogram of the yearly maximum mixed layer depth (identified as the first depth below which $|\sigma_\theta(z) - \sigma_\theta(z = -10 \text{ m})| > 0.01$, where $\sigma_\theta$ is the potential density referenced to sea level) in the deep water formation region, between $y = 45^\circ$ N and the northern boundary (i.e., the region north of $S_3$). We observe that $\Psi_{\text{MOC}}$ displays a structure consistent with previous works of the double-gyre configuration (Couespel et al., 2021; Ruan et al., 2023). Relative to the model truth R12, the R4 calculation has a rather large positive bias in the subpolar gyre region (Fig. 4b) and deep bias in the mixed layer (Fig. 4f), while the GEOM calculation has the converse (Fig. 4d, h). The observations are consistent with the expectation that the explicit eddy-mean interaction is too weak in R4, leading to a mixed layer that is too deep (since eddies are not able to counter the deepening of the mixed layer as much) and a MOC that is too strong. In GEOM the parameterised eddies lead to a response going too far the other way, leading to too shallow a mixed layer and too weak a MOC, reminiscent of the coarse resolution calculations reported in Fig. 3 and 4 of Ruan et al. (2023). However, we note that GEOM possesses a distribution in the mixed layer depth that is closer to R12, although there is a shallow bias.

The use of splitting appears to reduce the associated biases in the MOC and the distribution of sample isopycnals as seen in Fig. 4c, and reduce the deep biases of the mixed layer somewhat as seen in Fig. 4f, g (although the distribution is still rather wide compared to the distribution of the model truth in Fig. 4e). More quantitatively, the area-weighted average root-mean-square mismatch to R12 within the whole domain is of 1.21, 1.09, and 1.98 Sv in R4, SPLIT, and GEOM respectively. The median of the maximum mixed layer depth north of $S_3$ is 504, 765, 664 and 311 m, and the inter-quartile range is 232, 439, 401 and 197 m respectively for R12, R4, SPLIT, and GEOM.

One consequence of biases in the MOC is reflected in the ocean heat transports, diagnosed as

$$ \begin{aligned} \text{OHT} &= \rho_0 C_p \int_{z=-700}^{z=0} \int_{x=0}^{x=L_x} \overline{u\Theta} \, dx \, dz \\ &= \rho_0 C_p \int_{z=-700}^{z=0} \int_{x=0}^{x=L_x} \overline{(\bar{u} + u' + u^*)} (\overline{\Theta} + \Theta') \, dx \, dz, \end{aligned} \quad (9) $$

where $\rho_0$ is reference density at 1026 kg m$^{-3}$, $C_p = 3991.86 \text{ J K}^{-1}$ is the heat capacity, and $\Theta$ would be the Conservative Temperature (although we use a linear equation of state here). Fig. 5 shows the total meridional heat transport and vertical transport, zonally and vertically integrated over the top 700 m; we neglect the deeper parts of the ocean since these regions have not equilibrated within the 70 model year period we are considering. We note that the dominant contribution to the transports shown in Fig. 5 is from the mean component $\overline{u\Theta}$ (not shown). However, while the eddy components are subdominant in the overall transport, they are absolutely crucial for shaping the mean state and impacting $\overline{u\Theta}$.

The diagnosed meridional heat transport is mostly towards the north in the analysis region (within $S_1$ and $S_3$) and peaks near the Western Boundary Current (Fig. 5a). Relative to the model truth R12, R4 and GEOM possess a meridional heat transport that is too strong and too weak respectively, while SPLIT is much closer to model truth; quantitatively, the area-weighted average northward heat transport is 3.89, 7.29, 4.41, and -0.41 PW respectively for R12, R4, SPLIT and GEOM in the analysis domain.<fnref n="1" /> The improvements appear to come from a better representation of the stratification, which impacts both the transport and the heat content distribution.

For the vertical heat transport, there is a notable region with strong downward heat transport corresponding to the location of the Western Boundary Current (cf. Fig. 5b). The local biases between R12 with R4 and SPLIT are from the Western Boundary Current separating at a different latitude (there is almost no explicit representation of the Western Boundary Current in the GEOM case). The observation is partly the reason for the choice of a sample section $S_2$ at $y = 38^\circ$ N to be sufficiently away from the model Western Boundary Current, so that the associated section possesses a dynamical regime that is more comparable between calculations. Nevertheless, we can see that the magnitude of vertical heat transports of SPLIT is less than that of R4 over most of the region, and seems to be visually closer to the model truth R12 than GEOM. Quantitatively, the area-weighted average upward heat transport over the analysis region is 2.76, 5.74, 4.24 and $1.15 \times 10^{-3}$ PW respectively for R12, R4, SPLIT and GEOM, while the analogous downward heat transport is 3.44, 4.63, 3.39 and $0.84 \times 10^{-3}$ PW (the total heat transport is a small residual of the two, and is negative for R12). The key observation here is that the associated values for SPLIT are closer to the model truth R12, smaller than R4, and larger than that GEOM, consistent with our theoretical expectations.

<footnote n="1">¹ Note the values in Ruan et al. (2023) are smaller, because those are integrated over the whole domain depth. This was an inconsistent choice on our part.</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 3. The barotropic streamfunction $\Psi_{\text{baro}}$ (in Sv = $10^6 \text{ m}^3 \text{ s}^{-1}$) of (a) R12, (b) R4, (c) SPLIT, and (d) GEOM, with the zero contour added. The black dashed lines represent the sample latitudinal cross sections $S_{1,2,3}$ going from south to north.</figcaption></figure>

the subpolar gyre. While the works of Couespel et al. (2021) and Ruan et al. (2023) focus on the analysis within the subpolar gyre (defined over fixed geographical locations) primarily because it is the most bio-active region in the modelled domain, in this work we consider the whole domain but exclude the Fofonoff gyres as the analysis region when computing averaged/integrated quantities. For computing fluxes, we consider the cross sections at $y = 25^\circ$ and $45^\circ$ N marked on as black dashed lines in the figures as appropriate (denoted $S_{1,3}$), roughly as the southern boundary of the subtropical gyre, and the northern boundary of the subpolar gyre respectively. For completeness, we also mark on and compute fluxes over the $y = 38^\circ$ N section (denoted $S_2$), which is an empirically determined location that is sufficiently north of the simulated Western Boundary Current (see Fig. 4a-d). The analysis region of primary interest is the region bounded in the horizontal by $S_1$ and $S_3$, and in the vertical by $z = -700 \text{ m}$; we neglect the deeper parts of the ocean since these regions would presumably not have equilibrated within the 70 model year period we are considering.

As another measure of the circulation, we show in Fig. 4 the residual meridional overturning circulation (MOC) streamfunction $\Psi_{\text{MOC}}$, diagnosed as

$$ \Psi_{\text{MOC}} = \int_{z=-H}^{z=H} \int_{0}^{L_x} [\bar{v}(x, y, \bar{z}) + v^*(x, y, \bar{z})] \, dx \, d\bar{z}, \quad (8) $$

as well as a histogram of the yearly maximum mixed layer depth (identified as the first depth below which $|\sigma_\theta(z) - \sigma_\theta(z = -10 \text{ m})| > 0.01$, where $\sigma_\theta$ is the potential density referenced to sea level) in the deep water formation region, between $y = 45^\circ$ N and the northern boundary (i.e., the region north of $S_3$). We observe that $\Psi_{\text{MOC}}$ displays a structure consistent with previous works of the double-gyre configuration (Couespel et al., 2021; Ruan et al., 2023). Relative to the model truth R12, the R4 calculation has a rather large positive bias in the subpolar gyre region (Fig. 4b) and deep bias in the mixed layer (Fig. 4f), while the GEOM calculation has the converse (Fig. 4d, h). The observations are consistent with the expectation that the explicit eddy-mean interaction is too weak in R4, leading to a mixed layer that is too deep (since eddies are not able to counter the deepening of the mixed layer as much) and a MOC that is too strong. In GEOM the parameterised eddies lead to a response going too far the other way, leading to too shallow a mixed layer and too weak a MOC, reminiscent of the coarse resolution calculations reported in Fig. 3 and 4 of Ruan et al. (2023). However, we note that GEOM possesses a distribution in the mixed layer depth that is closer to R12, although there is a shallow bias.

The use of splitting appears to reduce the associated biases in the MOC and the distribution of sample isopycnals as seen in Fig. 4c, and reduce the deep biases of the mixed layer somewhat as seen in Fig. 4f, g (although the distribution is still rather wide compared to the distribution of the model truth in Fig. 4e). More quantitatively, the area-weighted average root-mean-square mismatch to R12 within the whole domain is of 1.21, 1.09, and 1.98 Sv in R4, SPLIT, and GEOM respectively. The median of the maximum mixed layer depth north of $S_3$ is 504, 765, 664 and 311 m, and the inter-quartile range is 232, 439, 401 and 197 m respectively for R12, R4, SPLIT, and GEOM.

One consequence of biases in the MOC is reflected in the ocean heat transports, diagnosed as

$$ \begin{aligned} \text{OHT} &= \rho_0 C_p \int_{z=-700}^{z=0} \int_{x=0}^{x=L_x} \overline{u\Theta} \, dx \, dz \\ &= \rho_0 C_p \int_{z=-700}^{z=0} \int_{x=0}^{x=L_x} \overline{(\bar{u} + u' + u^*)} (\overline{\Theta} + \Theta') \, dx \, dz, \end{aligned} \quad (9) $$

where $\rho_0$ is reference density at 1026 kg m$^{-3}$, $C_p = 3991.86 \text{ J K}^{-1}$ is the heat capacity, and $\Theta$ would be the Conservative Temperature (although we use a linear equation of state here). Fig. 5 shows the total meridional heat transport and vertical transport, zonally and vertically integrated over the top 700 m; we neglect the deeper parts of the ocean since these regions have not equilibrated within the 70 model year period we are considering. We note that the dominant contribution to the transports shown in Fig. 5 is from the mean component $\overline{u\Theta}$ (not shown). However, while the eddy components are subdominant in the overall transport, they are absolutely crucial for shaping the mean state and impacting $\overline{u\Theta}$.

The diagnosed meridional heat transport is mostly towards the north in the analysis region (within $S_1$ and $S_3$) and peaks near the Western Boundary Current (Fig. 5a). Relative to the model truth R12, R4 and GEOM possess a meridional heat transport that is too strong and too weak respectively, while SPLIT is much closer to model truth; quantitatively, the area-weighted average northward heat transport is 3.89, 7.29, 4.41, and -0.41 PW respectively for R12, R4, SPLIT and GEOM in the analysis domain.<fnref n="1" /> The improvements appear to come from a better representation of the stratification, which impacts both the transport and the heat content distribution.

For the vertical heat transport, there is a notable region with strong downward heat transport corresponding to the location of the Western Boundary Current (cf. Fig. 5b). The local biases between R12 with R4 and SPLIT are from the Western Boundary Current separating at a different latitude (there is almost no explicit representation of the Western Boundary Current in the GEOM case). The observation is partly the reason for the choice of a sample section $S_2$ at $y = 38^\circ$ N to be sufficiently away from the model Western Boundary Current, so that the associated section possesses a dynamical regime that is more comparable between calculations. Nevertheless, we can see that the magnitude of vertical heat transports of SPLIT is less than that of R4 over most of the region, and seems to be visually closer to the model truth R12 than GEOM. Quantitatively, the area-weighted average upward heat transport over the analysis region is 2.76, 5.74, 4.24 and $1.15 \times 10^{-3}$ PW respectively for R12, R4, SPLIT and GEOM, while the analogous downward heat transport is 3.44, 4.63, 3.39 and $0.84 \times 10^{-3}$ PW (the total heat transport is a small residual of the two, and is negative for R12). The key observation here is that the associated values for SPLIT are closer to the model truth R12, smaller than R4, and larger than that GEOM, consistent with our theoretical expectations.

<footnote n="1">¹ Note the values in Ruan et al. (2023) are smaller, because those are integrated over the whole domain depth. This was an inconsistent choice on our part.</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 4. (Top row) The Meridional Overturning Circulation streamfunction $\Psi_{MOC}$ (shading, in Sv = $10^6$ m$^3$ s$^{-1}$) and lines of constant potential density referenced to sea level (contours, in kg m$^{-3}$ with 0.5 kg m$^{-3}$ interval, for (a) R12, (b) R4, (c) SPLIT, and (d) GEOM. The black dashed lines represent the sample latitudinal cross sections $S_{1,2,3}$ going from south to north. (Bottom row) The histogram of yearly maximum mixed layer depth distributions (m, identified as the first depth below which $|\sigma_\theta(z) - \sigma_\theta(z = -10 \text{ m})| > 0.01$ where $\sigma_\theta$ is the potential density referenced to sea level, with 20 bins ranging from 0–1500 m) and median (indicated by horizontal blue line) over the northern area where deep water formation occurs, for (e) R12, (f) R4, (g) SPLIT, and (h) GEOM; the median and inter-quartile range (IQR) is shown (in units of m), and the axes of the histograms have been flipped for convenient visual comparison.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Total (mean and eddy) heat transport (in units of PW = $10^{15}$ W), integrated zonally and vertically over the top 700 m, for model truth (R12, black line) and eddy-permitting simulations (R4, blue line; SPLIT, red line; GEOM, yellow line). (a) Meridional heat transport $\overline{v\theta}$, positive values denoting northward transport. (b) Vertical heat transport $\overline{w\theta}$, positive values denoting upward transport. The black dashed lines represent the sample latitudinal cross sections $S_{1,2,3}$ going from south to north. We make a note that the values here are the transports integrated over the top 700 m, and are larger than the analogous values in Ruan et al. (2023), where the associated transports are integrated over the whole model depth; that former was an inconsistent choice on our part.</figcaption></figure>

## 3.2. Biogeochemical responses

Given the improvements to the physical responses in SPLIT achieved by employing the GEOMETRIC parameterisation with the splitting procedure, we might expect to observe related improvements in the biogeochemical responses. We focus on nitrate ($NO_3^-$), which contributes primarily to the NPP in the present set up (e.g., Couespel et al., 2021; Ruan et al., 2023). The improvement in the biogeochemical response can be seen from an improvement to the resulting nitrate distributions. The nitraclines, which largely mimics the isopycnal distribution since we expect transport to be constrained to along-isopycnal directions under the geostrophic assumption, is closer in SPLIT to the model truth (not shown for brevity, but see Fig. 4a-d for isopycnal distribution; cf. Ruan et al. 2023, Fig. 3a-c and Fig. 6a-c). The area-weighted average of $NO_3^-$ within the analysis domain are 12.02, 11.49, 11.81, and 10.55 mmol N m$^{-3}$ respectively for R12, R4, SPLIT and GEOM (where N is the nitrogen currency), demonstrating an improvement of SPLIT over R4 and GEOM in a way that is consistent with our expectations.

In Fig. 6 we show diagnostics related to the vertically integrated Net Primary Production (NPP; in units of mmol N m$^{-2}$ day$^{-1}$). Fig. 6a shows the distribution in the horizontal, and only the R12 simulation is shown since the general pattern of productive subpolar and oligotrophic subtropical gyre is similar in the eddy-permitting simulations. Fig. 6b shows the zonally averaged latitudinal distribution of the same quantity but over the set of calculations. Both R4 and SPLIT capture the general shape of the distribution for the R12 calculation, while

the GEOM calculation is too small in general, consistent with results from coarse resolution calculations (e.g., Fig. 5 of Ruan et al., 2023). Although the diagnosed NPP is weaker in SPLIT compared to R4 in the subpolar gyre region, the overshoot in the Western Boundary Current region is alleviated in SPLIT and is closer to the R12 model truth. Quantitatively, over the analysis region, the area-weighted average NPP is 2.64, 2.83, 2.61, and 1.31 mmol N m$^{-2}$ day$^{-1}$ respectively for R12, R4, SPLIT and GEOM. The diagnosed NPP for R4 and GEOM is too large and too small respectively, and SPLIT results in a much closer NPP to the model truth. The results in the biogeochemical response are consistent with our expected and diagnosed physical response, so there is evidence that we are getting an improved biogeochemical response because of a better physical response (cf. a case reported in Ruan et al. 2023, where one could obtain a reasonable biogeochemical response without necessarily having a consistent physical response).

Similar to the previous studies of Couespel et al. (2021) and Ruan et al. (2023), we diagnose the nitrate fluxes from advective and/or diffusive processes. The total nitrate advection is given by $\nabla \cdot (uN)$, where N denotes the nitrate concentration, and by the divergence theorem, the total supply in and out of the analysis region is expressed by

$$ \int_{\text{domain}} \nabla \cdot \overline{uN} \, dx \, dy \, dz = \left( \int_{S_1} + \int_{S_3} \right) \overline{vN} \, dx \, dz + \int_{A} \overline{wN} \, dx \, dy, \quad (10) $$

where A is the horizontal area between $S_1$ and $S_3$ at fixed height $z = z_0$, and we assume there is no surface input of N. The contributions can
<figure><img src="image_6.png" /><figcaption>Fig. 4. (Top row) The Meridional Overturning Circulation streamfunction $Ψ_{MOC}$ (shading, in Sv = $10^6$ m$^3$ s$^{-1}$) and lines of constant potential density referenced to sea level (contours, in kg m$^{-3}$ with 0.5 kg m$^{-3}$ interval, for (a) R12, (b) R4, (c) SPLIT, and (d) GEOM. The black dashed lines represent the sample latitudinal cross sections $S_{1,2,3}$ going from south to north. (Bottom row) The histogram of yearly maximum mixed layer depth distributions (m, identified as the first depth below which $|σ_θ(z) - σ_θ(z = -10 	ext{ m})| > 0.01$ where $σ_θ$ is the potential density referenced to sea level, with 20 bins ranging from 0–1500 m) and median (indicated by horizontal blue line) over the northern area where deep water formation occurs, for (e) R12, (f) R4, (g) SPLIT, and (h) GEOM; the median and inter-quartile range (IQR) is shown (in units of m), and the axes of the histograms have been flipped for convenient visual comparison.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Total (mean and eddy) heat transport (in units of PW = $10^{15}$ W), integrated zonally and vertically over the top 700 m, for model truth (R12, black line) and eddy-permitting simulations (R4, blue line; SPLIT, red line; GEOM, yellow line). (a) Meridional heat transport $\overline{v\theta}$, positive values denoting northward transport. (b) Vertical heat transport $\overline{w\theta}$, positive values denoting upward transport. The black dashed lines represent the sample latitudinal cross sections $S_{1,2,3}$ going from south to north. We make a note that the values here are the transports integrated over the top 700 m, and are larger than the analogous values in Ruan et al. (2023), where the associated transports are integrated over the whole model depth; that former was an inconsistent choice on our part.</figcaption></figure>

## 3.2. Biogeochemical responses

Given the improvements to the physical responses in SPLIT achieved by employing the GEOMETRIC parameterisation with the splitting procedure, we might expect to observe related improvements in the biogeochemical responses. We focus on nitrate ($NO_3^-$), which contributes primarily to the NPP in the present set up (e.g., Couespel et al., 2021; Ruan et al., 2023). The improvement in the biogeochemical response can be seen from an improvement to the resulting nitrate distributions. The nitraclines, which largely mimics the isopycnal distribution since we expect transport to be constrained to along-isopycnal directions under the geostrophic assumption, is closer in SPLIT to the model truth (not shown for brevity, but see Fig. 4a-d for isopycnal distribution; cf. Ruan et al. 2023, Fig. 3a-c and Fig. 6a-c). The area-weighted average of $NO_3^-$ within the analysis domain are 12.02, 11.49, 11.81, and 10.55 mmol N m$^{-3}$ respectively for R12, R4, SPLIT and GEOM (where N is the nitrogen currency), demonstrating an improvement of SPLIT over R4 and GEOM in a way that is consistent with our expectations.

In Fig. 6 we show diagnostics related to the vertically integrated Net Primary Production (NPP; in units of mmol N m$^{-2}$ day$^{-1}$). Fig. 6a shows the distribution in the horizontal, and only the R12 simulation is shown since the general pattern of productive subpolar and oligotrophic subtropical gyre is similar in the eddy-permitting simulations. Fig. 6b shows the zonally averaged latitudinal distribution of the same quantity but over the set of calculations. Both R4 and SPLIT capture the general shape of the distribution for the R12 calculation, while

the GEOM calculation is too small in general, consistent with results from coarse resolution calculations (e.g., Fig. 5 of Ruan et al., 2023). Although the diagnosed NPP is weaker in SPLIT compared to R4 in the subpolar gyre region, the overshoot in the Western Boundary Current region is alleviated in SPLIT and is closer to the R12 model truth. Quantitatively, over the analysis region, the area-weighted average NPP is 2.64, 2.83, 2.61, and 1.31 mmol N m$^{-2}$ day$^{-1}$ respectively for R12, R4, SPLIT and GEOM. The diagnosed NPP for R4 and GEOM is too large and too small respectively, and SPLIT results in a much closer NPP to the model truth. The results in the biogeochemical response are consistent with our expected and diagnosed physical response, so there is evidence that we are getting an improved biogeochemical response because of a better physical response (cf. a case reported in Ruan et al. 2023, where one could obtain a reasonable biogeochemical response without necessarily having a consistent physical response).

Similar to the previous studies of Couespel et al. (2021) and Ruan et al. (2023), we diagnose the nitrate fluxes from advective and/or diffusive processes. The total nitrate advection is given by $\nabla \cdot (uN)$, where N denotes the nitrate concentration, and by the divergence theorem, the total supply in and out of the analysis region is expressed by

$$ \int_{\text{domain}} \nabla \cdot \overline{uN} \, dx \, dy \, dz = \left( \int_{S_1} + \int_{S_3} \right) \overline{vN} \, dx \, dz + \int_{A} \overline{wN} \, dx \, dy, \quad (10) $$

where A is the horizontal area between $S_1$ and $S_3$ at fixed height $z = z_0$, and we assume there is no surface input of N. The contributions can
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Vertically integrated Net Primary Production (NPP, mmol N m⁻² day⁻¹, where N is the nitrogen currency). (a) Horizontal distribution for R12, and the distribution pattern in other simulations (R4, SPLIT, GEOM) is similar and thus is omitted. (b) The zonally averaged latitudinal distribution of vertically integrated NPP, for the model truth (R12, black line) and eddy-permitting simulations (R4, blue line; SPLIT, red line; GEOM, yellow line). The black dashed lines represent the sample latitudinal cross sections S<sub>1,2,3</sub> going from south to north.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Total advective nitrate supply (mmol N day⁻¹, where N is the nitrogen currency unit) for the model truth (R12, black line) and eddy-permitting simulations (R4, blue line; SPLIT, red line; GEOM, yellow line). (a) The vertically varying vertical advective supply horizontally integrated. (b, c, d) The meridional advective contribution as a cumulative vertical integral at sample sections S<sub>1,2,3</sub>. The contribution at S<sub>3</sub> is calculated with an extra minus sign, so negative values in panel d indicate a flux out of the analysis region. Also note that in panel c the blue line overlaps with the red line.</figcaption></figure>

further be decomposed into explicit and parameterised eddy compo-
nents as in Eq. (9). In Fig. 7 we show the total nitrate fluxes, where the
vertical flux is diagnosed by performing a horizontal integral (panel
*a*), while the meridional fluxes across S<sub>1,2,3</sub> are the zonal integral over
S<sub>1,2,3</sub> that are then cumulatively integrated in the vertical (panels b, *c*, *d*).
We note that while the total is shown, the dominant contribution is from
the mean component **v**̅*N* and **w**̅*N*, rather than from the eddy component
(explicit and/or parameterised), consistent with the previous results
reported in Couespel et al. (2021) and Ruan et al. (2023). Vertical
diffusion is large over the top 50 m or so, while lateral diffusion is of
secondary importance over all depths (not shown). Again, while the
eddy contribution to tracer transport may be of secondary importance,
the eddies are crucial in shaping the mean stratification, which ends
up dictating the overall large-scale supply of nitrate.

The nitrate advection profiles of R4 and SPLIT agree reasonably
well with R12, with SPLIT being generally of smaller magnitude than
R4, and certainly an improvement on GEOM (and analogous coarse
resolution calculations, such as those in Fig. 7 of Ruan et al. 2023). The
main difference appears in the vertical nitrate supply in Fig. 7a, where
the vertical supply of R4 and SPLIT are too small near the upper parts of
the modelled ocean (and, interestingly, the diagnosed values of GEOM
agree better with R12 here), while they are larger at the deeper regions.
The smaller supply of nitrate particularly in the vertical, is presumably
the dominant contribution to why the NPP in SPLIT is smaller than R4,
and closer to that of the model truth R12.

A final comment we make is that there are large differences in the
associated meridional transport if S₂ is within the explicitly represented
Western Boundary Current. For comparison reasons, we chose S₂ to
be north of all the explicitly represented Western Boundary Currents
across the set of simulations, and somewhere near the southern bound-
ary of the modelled subpolar gyre. While it may be possible to compute
the flux across some contour related to the Western Boundary Current
defined dynamically, for simplicity reasons we have opted to simply
choose a sample section to provide a representative diagnostic.

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Raw difference of the $\kappa_{gm}(x, y)$ (m²s⁻¹) distribution between the climate change scenario and the control scenario (see Fig. 2) for (a) SPLIT, and (b) GEOM. The black dashed lines represent the sample latitudinal cross sections $S_{1,2,3}$ going from south to north.</figcaption></figure>

4. Noteworthy characteristics of the sensitivities under idealised climate change

4.1. *Physical responses*

In the previous work of Ruan et al. (2023) we investigated the
performance of GEOMETRIC under idealised climate change scenarios
for coarse resolution models, and found that the use of GEOMETRIC
improved on the sensitivities, at least compared to the standard GM
scheme with a constant $\kappa_{gm}$. We perform a similar set of calculations un-
der the same idealised climate change scenario detailed in Section 2.4
to investigate the responses between the eddy-permitting models R4,
SPLIT and GEOM.

Fig. 8 shows the raw difference of κ<sub>gm</sub>(x, y) between the climate
change and control case (cf. Fig. 2) for SPLIT and GEOM. The change
in κ<sub>gm</sub> mainly appears near the Western Boundary Current region and
the northern boundary, and is relatively modest for SPLIT compared to
GEOM, because it is the large-scale filtered stratification that is used in
the calculation of the parameterised total eddy energy (i.e., Eq. (5) and
(6)). The change in κ<sub>gm</sub> in GEOM is large in the deep water forming
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Vertically integrated Net Primary Production (NPP, mmol N m⁻² day⁻¹, where N is the nitrogen currency). (a) Horizontal distribution for R12, and the distribution pattern in other simulations (R4, SPLIT, GEOM) is similar and thus is omitted. (b) The zonally averaged latitudinal distribution of vertically integrated NPP, for the model truth (R12, black line) and eddy-permitting simulations (R4, blue line; SPLIT, red line; GEOM, yellow line). The black dashed lines represent the sample latitudinal cross sections S<sub>1,2,3</sub> going from south to north.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Total advective nitrate supply (mmol N day⁻¹, where N is the nitrogen currency unit) for the model truth (R12, black line) and eddy-permitting simulations (R4, blue line; SPLIT, red line; GEOM, yellow line). (a) The vertically varying vertical advective supply horizontally integrated. (b, c, d) The meridional advective contribution as a cumulative vertical integral at sample sections S<sub>1,2,3</sub>. The contribution at S<sub>3</sub> is calculated with an extra minus sign, so negative values in panel d indicate a flux out of the analysis region. Also note that in panel c the blue line overlaps with the red line.</figcaption></figure>

further be decomposed into explicit and parameterised eddy compo-
nents as in Eq. (9). In Fig. 7 we show the total nitrate fluxes, where the
vertical flux is diagnosed by performing a horizontal integral (panel
*a*), while the meridional fluxes across S<sub>1,2,3</sub> are the zonal integral over
S<sub>1,2,3</sub> that are then cumulatively integrated in the vertical (panels b, *c*, *d*).
We note that while the total is shown, the dominant contribution is from
the mean component *v*̅N and *w*̅N, rather than from the eddy component
(explicit and/or parameterised), consistent with the previous results
reported in Couespel et al. (2021) and Ruan et al. (2023). Vertical
diffusion is large over the top 50 m or so, while lateral diffusion is of
secondary importance over all depths (not shown). Again, while the
eddy contribution to tracer transport may be of secondary importance,
the eddies are crucial in shaping the mean stratification, which ends
up dictating the overall large-scale supply of nitrate.

The nitrate advection profiles of R4 and SPLIT agree reasonably well with R12, with SPLIT being generally of smaller magnitude than R4, and certainly an improvement on GEOM (and analogous coarse resolution calculations, such as those in Fig. 7 of Ruan et al. 2023). The main difference appears in the vertical nitrate supply in Fig. 7a, where the vertical supply of R4 and SPLIT are too small near the upper parts of the modelled ocean (and, interestingly, the diagnosed values of GEOM agree better with R12 here), while they are larger at the deeper regions. The smaller supply of nitrate particularly in the vertical, is presumably the dominant contribution to why the NPP in SPLIT is smaller than R4, and closer to that of the model truth R12.

A final comment we make is that there are large differences in the
associated meridional transport if S₂ is within the explicitly represented
Western Boundary Current. For comparison reasons, we chose S₂ to
be north of all the explicitly represented Western Boundary Currents
across the set of simulations, and somewhere near the southern bound-
ary of the modelled subpolar gyre. While it may be possible to compute
the flux across some contour related to the Western Boundary Current
defined dynamically, for simplicity reasons we have opted to simply
choose a sample section to provide a representative diagnostic.

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Raw difference of the $\kappa_{gm}(x, y)$ (m²s⁻¹) distribution between the climate change scenario and the control scenario (see Fig. 2) for (a) SPLIT, and (b) GEOM. The black dashed lines represent the sample latitudinal cross sections $S_{1,2,3}$ going from south to north.</figcaption></figure>

4. Noteworthy characteristics of the sensitivities under idealised climate change

**4.1. Physical responses**

In the previous work of Ruan et al. (2023) we investigated the
performance of GEOMETRIC under idealised climate change scenarios
for coarse resolution models, and found that the use of GEOMETRIC
improved on the sensitivities, at least compared to the standard GM
scheme with a constant $\kappa_{gm}$. We perform a similar set of calculations un-
der the same idealised climate change scenario detailed in Section 2.4
to investigate the responses between the eddy-permitting models R4,
SPLIT and GEOM.

Fig. 8 shows the raw difference of κ<sub>gm</sub>(x, y) between the climate
change and control case (cf. Fig. 2) for SPLIT and GEOM. The change
in κ<sub>gm</sub> mainly appears near the Western Boundary Current region and
the northern boundary, and is relatively modest for SPLIT compared to
GEOM, because it is the large-scale filtered stratification that is used in
the calculation of the parameterised total eddy energy (i.e., Eq. (5) and
(6)). The change in κ<sub>gm</sub> in GEOM is large in the deep water forming
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** (Top row) Raw difference of the meridional overturning streamfunction $\Psi_{\text{MOC}}$ (Sv = $10^6$ m$^3$ s$^{-1}$) between the climate change and the control scenario (Fig. 4a-d), for (a) R12, (b) R4, (c) SPLIT, and (d) GEOM; negative values largely indicate a weakening of the overturning strength. The black dashed lines represent the sample latitudinal cross sections $S_{1,2,3}$ going from south to north. (Bottom row) The histogram of yearly maximum mixed layer depth distributions ($m$, identified as the first depth below which $|\sigma_{\theta}(z) - \sigma_{\theta}(z = -10 \text{ m})| > 0.01$ where $\sigma_{\theta}$ is the potential density referenced to sea level, with 20 bins ranging from 0–1500 m) and median (indicated by the horizontal line) over the area where deep water formation occurs, for the climate change scenario (in red) and control scenario (in blue, see Fig. 4e-h), for (e) R12, (f) R4, (g) SPLIT, and (h) GEOM; the median and inter-quartile range (IQR) is shown (in units of $m$), and the axes of the histograms have been flipped for convenient visual comparison.</figcaption></figure>

region towards the northern boundary, consistent with but larger in magnitude than the corresponding result in Fig. 8b of Ruan et al. (2023) for the coarse resolution calculation. As reported previously in Ruan et al. (2023), this significant increase in $\kappa_{\text{gm}}$ in the deep water formation region affects the mixed layer depth and extent of deep water formation, which affects the overall stratification in the domain, impacting the meridional overturning circulation.

Fig. 9a-d shows the raw difference between the meridional overturning streamfunction $\Psi_{\text{MOC}}$ under the climate change and control scenario, while Fig. 9e-h shows the histogram of the yearly maximum mixed layer depth distributions under both scenarios (cf. Fig. 4). Overall, the overturning circulation weakens and the mixed layer depth shoals, consistent with the results from the previous works of Couespel et al. (2021) and Ruan et al. (2023). Note that the change of $\Psi_{\text{MOC}}$ in R12 is somewhat different to that reported in Fig. 9a of Ruan et al. (2023): the strength in the subtropical gyre here *decreases*, while it *increases* in the previous works, presumably related to the presence of the non-trivial bathymetry in this work. The alternating positive-negative pattern near the centre of the domain corresponds to a shift of the Western Boundary Current northwards (corresponding to a shift in the purple pattern in Fig. 4a). We note that the shift in the Western Boundary Current as seen in $\Psi_{\text{MOC}}$ is more noticeable in SPLIT. The decrease in the subpolar gyre is less in SPLIT relative to R4 and closer to R12 visually, although R4 possesses a change in the subtropical gyre that is closer to R12. R4 and SPLIT still have a slightly deep bias in the mixed layer, but overall the shift is not unreasonable compared to R12. By contrast, GEOM has a weak change in the $\Psi_{\text{MOC}}$, but only because the control $\Psi_{\text{MOC}}$ in Fig. 4d is already rather small. Additionally, the mixed layer in GEOM is too shallow, as expected from the strong increase in $\kappa_{\text{gm}}$ in the deep water forming region, consistent with results from Ruan et al. (2023).

The changes in both the meridional and vertical heat transport are consistent with the observed changes in $\Psi_{\text{MOC}}$, and are perhaps best quantified by simply stating the diagnosed values averaged with the analysis region. The area-weighted averaged northward heat transport is 4.67, 5.05, 3.17 and 0.56 PW for R12, R4, SPLIT and GEOM respectively. In this setting, the diagnosed northward heat transport for R4 is now closer to R12 than that of SPLIT. Similarly, the total upward heat transport is 3.41, 6.46, 5.40 and $1.07 \times 10^{-3}$ PW for R12, R4, SPLIT and GEOM respectively, while for downward heat transport is 3.18, 4.96, 4.35 and $1.12 \times 10^{-3}$ PW; the values for SPLIT are more consistent with R12 and smaller in magnitude than R4 as expected. Note however the sensitivities of the meridional and vertical heat

transports of both R4 and SPLIT between climate change and control scenarios differ in magnitude and sometimes in sign (see Table 2), suggesting improvements to the sensitivities relative to the R12 model truth requires further investigation.

## 4.2. Biogeochemical responses

In Fig. 10 we show diagnostics in relation to the vertically integrated NPP. Fig. 10a shows the horizontal distribution of the raw difference between climate change and control scenario for R12, where we see there is an overall decline of NPP, in line with results from Couespel et al. (2021) and Ruan et al. (2023). Minor differences in the Western Boundary Current region arise compared to previous works, presumably because of the presence of the non-trivial bathymetry leading to a different representation of the circulation (e.g. Fig. 11a of Ruan et al. 2023). The corresponding figures for the eddy-permitting calculations have been omitted because the changes are qualitatively similar. Fig. 10b shows the analogous zonally averaged diagnostic across the set of experiments. We see that the decrease near the Western Boundary Current in SPLIT matches better with R12, and both SPLIT and R4 perform similarly in the subpolar gyre, both with a positive bias. The GEOM calculation displays substantial sensitivity in the subpolar gyre, because of the significant change in $\kappa_{\text{gm}}$ (Fig. 8b), analogous to the sensitivity reported in Ruan et al. (2023). The overall sensitivity within the analysis domain between R12, R4 and SPLIT are not significantly different. Quantitatively, the area-weighted vertically integrated NPP is 2.18, 2.36, 2.18 and $1.07 \text{ mmol N m}^{-2} \text{ day}^{-1}$, with a corresponding percentage difference of -17.2%, -16.5%, -16.6%, and -18.3% for R12, R4, SPLIT and GEOM respectively. While the SPLIT experiment arguably displays a better agreement with R12, the R4 experiment performs reasonably (and certainly better than GEOM in the present eddy-permitting regime).

In Fig. 11 we show the raw difference between the vertical advective nitrate supply and the meridional advective nitrate supply over the vertical planes $S_{1,2,3}$ (cf. Fig. 7). In general, the vertical nitrate supply and meridional nitrate supply in $S_{1,2}$ decreases, but increases in $S_3$ (which because of the minus sign in the calculation corresponds to a *negative* supply or positive flux out of the analysis domain), related to a weakening of the overturning circulation (cf. Fig. 9a-d). There is however a general disagreement in terms of the changes in the advective contributions in the eddy-permitting calculations with R12, which is in contrast to the coarse resolution results reported in Fig. 13 of Ruan et al. (2023). This could be due to the presence of bathymetry
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** (Top row) Raw difference of the meridional overturning streamfunction $\Psi_{\text{MOC}}$ (Sv = $10^6$ m$^3$ s$^{-1}$) between the climate change and the control scenario (Fig. 4a-d), for (a) R12, (b) R4, (c) SPLIT, and (d) GEOM; negative values largely indicate a weakening of the overturning strength. The black dashed lines represent the sample latitudinal cross sections $S_{1,2,3}$ going from south to north. (Bottom row) The histogram of yearly maximum mixed layer depth distributions ($m$, identified as the first depth below which $|\sigma_{\theta}(z) - \sigma_{\theta}(z = -10 \text{ m})| > 0.01$ where $\sigma_{\theta}$ is the potential density referenced to sea level, with 20 bins ranging from 0–1500 m) and median (indicated by the horizontal line) over the area where deep water formation occurs, for the climate change scenario (in red) and control scenario (in blue, see Fig. 4e-h), for (e) R12, (f) R4, (g) SPLIT, and (h) GEOM; the median and inter-quartile range (IQR) is shown (in units of $m$), and the axes of the histograms have been flipped for convenient visual comparison.</figcaption></figure>

region towards the northern boundary, consistent with but larger in magnitude than the corresponding result in Fig. 8b of Ruan et al. (2023) for the coarse resolution calculation. As reported previously in Ruan et al. (2023), this significant increase in $\kappa_{\text{gm}}$ in the deep water formation region affects the mixed layer depth and extent of deep water formation, which affects the overall stratification in the domain, impacting the meridional overturning circulation.

Fig. 9a-d shows the raw difference between the meridional overturning streamfunction $\Psi_{\text{MOC}}$ under the climate change and control scenario, while Fig. 9e-h shows the histogram of the yearly maximum mixed layer depth distributions under both scenarios (cf. Fig. 4). Overall, the overturning circulation weakens and the mixed layer depth shoals, consistent with the results from the previous works of Couespel et al. (2021) and Ruan et al. (2023). Note that the change of $\Psi_{\text{MOC}}$ in R12 is somewhat different to that reported in Fig. 9a of Ruan et al. (2023): the strength in the subtropical gyre here *decreases*, while it *increases* in the previous works, presumably related to the presence of the non-trivial bathymetry in this work. The alternating positive-negative pattern near the centre of the domain corresponds to a shift of the Western Boundary Current northwards (corresponding to a shift in the purple pattern in Fig. 4a). We note that the shift in the Western Boundary Current as seen in $\Psi_{\text{MOC}}$ is more noticeable in SPLIT. The decrease in the subpolar gyre is less in SPLIT relative to R4 and closer to R12 visually, although R4 possesses a change in the subtropical gyre that is closer to R12. R4 and SPLIT still have a slightly deep bias in the mixed layer, but overall the shift is not unreasonable compared to R12. By contrast, GEOM has a weak change in the $\Psi_{\text{MOC}}$, but only because the control $\Psi_{\text{MOC}}$ in Fig. 4d is already rather small. Additionally, the mixed layer in GEOM is too shallow, as expected from the strong increase in $\kappa_{\text{gm}}$ in the deep water forming region, consistent with results from Ruan et al. (2023).

The changes in both the meridional and vertical heat transport are consistent with the observed changes in $\Psi_{\text{MOC}}$, and are perhaps best quantified by simply stating the diagnosed values averaged with the analysis region. The area-weighted averaged northward heat transport is 4.67, 5.05, 3.17 and 0.56 PW for R12, R4, SPLIT and GEOM respectively. In this setting, the diagnosed northward heat transport for R4 is now closer to R12 than that of SPLIT. Similarly, the total upward heat transport is 3.41, 6.46, 5.40 and $1.07 \times 10^{-3}$ PW for R12, R4, SPLIT and GEOM respectively, while for downward heat transport is 3.18, 4.96, 4.35 and $1.12 \times 10^{-3}$ PW; the values for SPLIT are more consistent with R12 and smaller in magnitude than R4 as expected. Note however the sensitivities of the meridional and vertical heat

transports of both R4 and SPLIT between climate change and control scenarios differ in magnitude and sometimes in sign (see Table 2), suggesting improvements to the sensitivities relative to the R12 model truth requires further investigation.

## 4.2. Biogeochemical responses

In Fig. 10 we show diagnostics in relation to the vertically integrated NPP. Fig. 10a shows the horizontal distribution of the raw difference between climate change and control scenario for R12, where we see there is an overall decline of NPP, in line with results from Couespel et al. (2021) and Ruan et al. (2023). Minor differences in the Western Boundary Current region arise compared to previous works, presumably because of the presence of the non-trivial bathymetry leading to a different representation of the circulation (e.g. Fig. 11a of Ruan et al. 2023). The corresponding figures for the eddy-permitting calculations have been omitted because the changes are qualitatively similar. Fig. 10b shows the analogous zonally averaged diagnostic across the set of experiments. We see that the decrease near the Western Boundary Current in SPLIT matches better with R12, and both SPLIT and R4 perform similarly in the subpolar gyre, both with a positive bias. The GEOM calculation displays substantial sensitivity in the subpolar gyre, because of the significant change in $\kappa_{\text{gm}}$ (Fig. 8b), analogous to the sensitivity reported in Ruan et al. (2023). The overall sensitivity within the analysis domain between R12, R4 and SPLIT are not significantly different. Quantitatively, the area-weighted vertically integrated NPP is 2.18, 2.36, 2.18 and $1.07 \text{ mmol N m}^{-2} \text{ day}^{-1}$, with a corresponding percentage difference of -17.2%, -16.5%, -16.6%, and -18.3% for R12, R4, SPLIT and GEOM respectively. While the SPLIT experiment arguably displays a better agreement with R12, the R4 experiment performs reasonably (and certainly better than GEOM in the present eddy-permitting regime).

In Fig. 11 we show the raw difference between the vertical advective nitrate supply and the meridional advective nitrate supply over the vertical planes $S_{1,2,3}$ (cf. Fig. 7). In general, the vertical nitrate supply and meridional nitrate supply in $S_{1,2}$ decreases, but increases in $S_3$ (which because of the minus sign in the calculation corresponds to a *negative* supply or positive flux out of the analysis domain), related to a weakening of the overturning circulation (cf. Fig. 9a-d). There is however a general disagreement in terms of the changes in the advective contributions in the eddy-permitting calculations with R12, which is in contrast to the coarse resolution results reported in Fig. 13 of Ruan et al. (2023). This could be due to the presence of bathymetry
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Raw difference of the vertically integrated Net Primary Production (NPP, mmol N m⁻² day⁻¹, where N is the nitrogen currency) between the climate change scenario and the corresponding control scenario (see Fig. 6a). (a) Horizontal distribution for R12; the distribution pattern is similar for the eddy-permitting simulations and have been omitted. (b) The zonally averaged latitudinal distribution for model truth (R12, black line) and eddy-permitting simulations (R4, blue line; SPLIT, red line; GEOM, yellow line). The black dashed lines represent the sample latitudinal cross sections S<sub>1,2,3</sub> going from south to north.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Raw difference of the total advective nitrate supply (mmol N day⁻¹, where N is the nitrogen currency unit) into the analysis domain between the climate change scenario and the corresponding control scenario (see Fig. 7) for model truth (R12, black line) and eddy-permitting simulations (R4, blue line; SPLIT, red line; GEOM, yellow line). (a) The vertically varying vertical advective supply horizontally integrated. (b–d) The meridional advective contribution as a cumulative vertical integral at sample sections S<sub>1,2,3</sub>. The contribution at S<sub>3</sub> is calculated with an extra minus sign, so that positive values here indicate a positive supply (or decrease in the flux out of the analysis domain).</figcaption></figure>

in the present model or the absence of isoneutral diffusion, but also
could be that the statement “*a GM-based parameterisation turns an eddy-
permitting model into an expensive coarse resolution model*” as we have
hypothesised is overly simplified. While the overall integrated response
of the SPLIT calculation appears to be reasonable relative to the model
truth, further investigations and/or proposal for improvements are
warranted in due course.

5. Conclusions and discussions

With increasing computational power, eddy-permitting ocean mod-
els that can partially resolve mesoscale eddies become increasingly
feasible, and are tractable targets for the next generation of Earth
System Models. While it is known that such models offer some benefits
in reducing biases (e.g., Fox-Kemper et al., 2019; Hewitt et al., 2020;
Beech et al., 2022), it is also known that the eddy-mean interaction
are somewhat misrepresented, and some degree of mesoscale eddy
parameterisation may still be beneficial. Our work here contributes to
the examination of the effect of different ways of representing the eddy-
mean interaction on the associated model physical and biogeochemical
response in an idealised double-gyre model, based on the model and
experiment set up from the related previous works (Couespel et al.,
2021; Ruan et al., 2023). A difference in the model set up compared
with the previous works is that we employ a non-trivial bathymetry,
with a slope on the west and east sides (Jackson et al. 2006; cf. Eq. (4)
and Fig. 1a). The qualitative conclusions we draw from this work are
however robust also in the absence of model bathymetry (not shown).

We focus our attention on examining a new way of employing
a Gent–McWilliams (GM) scheme (Gent and McWilliams, 1990; Gent
et al., 1995) in eddy-permitting ocean models, termed splitting in the
previous work of Mak et al. (2023), to the present double-gyre system.
Focusing on the GM-based version GEOMETRIC (Marshall et al., 2012;
Mak et al., 2022b), we compare model responses of a model employing
GEOMETRIC with splitting against model responses in an eddy-rich
model truth, and eddy-permitting cases where no GM-based parameter-
isation is employed, and where a GM-based scheme is used as-is. We

expect that in the case with no parameterisation, the modelled stratification will be too deep and the meridional overturning circulation too strong, because the represented eddy-mean interaction in an eddy-permitting calculation is too weak. As a consequence, the modelled heat transports and biogeochemical response in the net primary production (NPP) will be too large relative to the model truth. On the other hand, a model using the GM-based version of GEOMETRIC as-is with no splitting is expected to largely behave like a coarse resolution model. We thus expect that the corresponding modelled meridional overturning circulation will be too weak, because the parameterisation used as-is ends up taking over and removing contributions from explicit eddies permitted by the model resolution (cf. Fig. 1 of Mak et al., 2023), and the associated modelled heat transport and NPP is too small. Our hypothesis is that GEOMETRIC with splitting is able to better capture physical and biogeochemical responses displayed by an eddy-rich model truth, and in a way that is physically consistent. The use of splitting has been shown to allow the explicit and parameterised eddy components to co-exist (Mak et al., 2023). The resulting responses are largely like that of an eddy-permitting calculation without parameterisation, but with an extra contribution affecting the mean state from the parameterisation, so that the modelled stratification, the meridional overturning circulation, and the magnitude of the modelled heat transports and NPP are closer to the model truth.

Table 2 summarises the metrics of interest in this work. The general
conclusion is that, indeed, the use of GEOMETRIC and splitting broadly
improves upon the modelled biases relative to the eddy-rich model
truth over both the eddy-permitting models with no parameterisation
or with the GM-based version of GEOMETRIC applied as-is, in the ex-
pected way detailed in the previous paragraph, under both the control
and idealised climate change scenario. We reiterate that, while the
more detailed analyses performed suggest that the eddy contribution
(explicit and/or parameterised) to tracer transport is rather small, the
eddies are essential for shaping the mean-state, which ultimately leads
to substantial changes in the bulk diagnostics (Couespel et al., 2021;
Ruan et al., 2023). It should be noted that the model sensitivities
could still be improved upon. While the use of splitting certainly
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Raw difference of the vertically integrated Net Primary Production (NPP, mmol N m⁻² day⁻¹, where N is the nitrogen currency) between the climate change scenario and the corresponding control scenario (see Fig. 6a). (a) Horizontal distribution for R12; the distribution pattern is similar for the eddy-permitting simulations and have been omitted. (b) The zonally averaged latitudinal distribution for model truth (R12, black line) and eddy-permitting simulations (R4, blue line; SPLIT, red line; GEOM, yellow line). The black dashed lines represent the sample latitudinal cross sections $S_{1,2,3}$ going from south to north.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Raw difference of the total advective nitrate supply (mmol N day⁻¹, where N is the nitrogen currency unit) into the analysis domain between the climate change scenario and the corresponding control scenario (see Fig. 7) for model truth (R12, black line) and eddy-permitting simulations (R4, blue line; SPLIT, red line; GEOM, yellow line). (a) The vertically varying vertical advective supply horizontally integrated. (b–d) The meridional advective contribution as a cumulative vertical integral at sample sections $S_{1,2,3}$. The contribution at $S_3$ is calculated with an extra minus sign, so that positive values here indicate a positive supply (or decrease in the flux out of the analysis domain).</figcaption></figure>

in the present model or the absence of isoneutral diffusion, but also could be that the statement “*a GM-based parameterisation turns an eddy-permitting model into an expensive coarse resolution model*” as we have hypothesised is overly simplified. While the overall integrated response of the SPLIT calculation appears to be reasonable relative to the model truth, further investigations and/or proposal for improvements are warranted in due course.

## 5. Conclusions and discussions

With increasing computational power, eddy-permitting ocean models that can partially resolve mesoscale eddies become increasingly feasible, and are tractable targets for the next generation of Earth System Models. While it is known that such models offer some benefits in reducing biases (e.g., Fox-Kemper et al., 2019; Hewitt et al., 2020; Beech et al., 2022), it is also known that the eddy-mean interaction are somewhat misrepresented, and some degree of mesoscale eddy parameterisation may still be beneficial. Our work here contributes to the examination of the effect of different ways of representing the eddy-mean interaction on the associated model physical and biogeochemical response in an idealised double-gyre model, based on the model and experiment set up from the related previous works (Couespel et al., 2021; Ruan et al., 2023). A difference in the model set up compared with the previous works is that we employ a non-trivial bathymetry, with a slope on the west and east sides (Jackson et al. 2006; cf. Eq. (4) and Fig. 1a). The qualitative conclusions we draw from this work are however robust also in the absence of model bathymetry (not shown).

We focus our attention on examining a new way of employing a Gent-McWilliams (GM) scheme (Gent and McWilliams, 1990; Gent et al., 1995) in eddy-permitting ocean models, termed splitting in the previous work of Mak et al. (2023), to the present double-gyre system. Focusing on the GM-based version GEOMETRIC (Marshall et al., 2012; Mak et al., 2022b), we compare model responses of a model employing GEOMETRIC with splitting against model responses in an eddy-rich model truth, and eddy-permitting cases where no GM-based parameterisation is employed, and where a GM-based scheme is used as-is. We

expect that in the case with no parameterisation, the modelled stratification will be too deep and the meridional overturning circulation too strong, because the represented eddy-mean interaction in an eddy-permitting calculation is too weak. As a consequence, the modelled heat transports and biogeochemical response in the net primary production (NPP) will be too large relative to the model truth. On the other hand, a model using the GM-based version of GEOMETRIC as-is with no splitting is expected to largely behave like a coarse resolution model. We thus expect that the corresponding modelled meridional overturning circulation will be too weak, because the parameterisation used as-is ends up taking over and removing contributions from explicit eddies permitted by the model resolution (cf. Fig. 1 of Mak et al., 2023), and the associated modelled heat transport and NPP is too small. Our hypothesis is that GEOMETRIC with splitting is able to better capture physical and biogeochemical responses displayed by an eddy-rich model truth, and in a way that is physically consistent. The use of splitting has been shown to allow the explicit and parameterised eddy components to co-exist (Mak et al., 2023). The resulting responses are largely like that of an eddy-permitting calculation without parameterisation, but with an extra contribution affecting the mean state from the parameterisation, so that the modelled stratification, the meridional overturning circulation, and the magnitude of the modelled heat transports and NPP are closer to the model truth.

Table 2 summarises the metrics of interest in this work. The general conclusion is that, indeed, the use of GEOMETRIC and splitting broadly improves upon the modelled biases relative to the eddy-rich model truth over both the eddy-permitting models with no parameterisation or with the GM-based version of GEOMETRIC applied as-is, in the expected way detailed in the previous paragraph, under both the control and idealised climate change scenario. We reiterate that, while the more detailed analyses performed suggest that the eddy contribution (explicit and/or parameterised) to tracer transport is rather small, the eddies are essential for shaping the mean-state, which ultimately leads to substantial changes in the bulk diagnostics (Couespel et al., 2021; Ruan et al., 2023). It should be noted that the model sensitivities could still be improved upon. While the use of splitting certainly
<caption>**Table 2**</caption>

Summary of diagnostics and their sensitivities for the set of calculations, all analysed within the analysis domain between $y = 25$ and $45^{\circ}$ N ($S_1$ and $S_3$ in the text) and $z = -700$ m, except for mixed layer depth, which is analysed north of $y = 45^{\circ}$ N in the deep water formation region. The bracketed numbers denote the percentage differences of the relevant diagnostic between the climate change and control scenario. $L^2$ sensitivity denotes the area-weighted average root-mean-square difference between the climate change and control scenario. A dash is given if there is no obvious evidence that SPLIT leads to an improvement over R4. Note that the values of the northward heat transport here are larger than those reported in Ruan et al. (2023), which was computed from an integral over the whole model depth, and was an inconsistent choice on our part.

<table><thead><tr><th>Diagnostic</th><th>R12 values</th><th>R4 values</th><th>SPLIT values</th><th>GEOM values</th><th>Improve by SPLIT</th></tr></thead><tbody><tr><td><strong>overturning circulation (Sv)</strong><br>(Figs. 4a-d and 9a-d)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>L<sup>2</sup> mismatch rel. R12 (CTL)</td><td>-</td><td>1.21</td><td>1.09</td><td>1.98</td><td>✓</td></tr><tr><td>L<sup>2</sup> mismatch rel. R12 (CC)</td><td>-</td><td>0.70</td><td>0.69</td><td>1.58</td><td>-</td></tr><tr><td><strong>northern mixed layer depth (m)</strong><br>(Figs. 4e-h and 9e-h)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>median (CTL)</td><td>504</td><td>765</td><td>664</td><td>311</td><td>✓</td></tr><tr><td>median (CC)</td><td>347 (-31.2%)</td><td>478 (-37.5%)</td><td>422 (-36.5%)</td><td>225 (-27.6%)</td><td>✓(✓)</td></tr><tr><td>quartile range (CTL)</td><td>232</td><td>439</td><td>401</td><td>197</td><td>✓</td></tr><tr><td>quartile range (CC)</td><td>111 (-52.4%)</td><td>199 (-54.6%)</td><td>172 (-57.1%)</td><td>128 (-34.8%)</td><td>✓(x)</td></tr><tr><td><strong>northward heat transport (PW)</strong><br>(Fig. 5a)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>area average (CTL)</td><td>3.89</td><td>7.29</td><td>4.41</td><td>-0.41</td><td>✓</td></tr><tr><td>area average (CC)</td><td>4.67 (+20.0%)</td><td>5.05 (-30.7%)</td><td>3.17 (-28.2%)</td><td>0.56 (-235.8%)</td><td>✓(-)</td></tr><tr><td>sensitivity (L<sup>2</sup>)</td><td>2.47</td><td>3.03</td><td>2.10</td><td>1.41</td><td>✓</td></tr><tr><td><strong>vertical heat transport (10<sup>-3</sup> PW)</strong><br>(Fig. 5b)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>area average upward (CTL)</td><td>2.76</td><td>5.74</td><td>4.24</td><td>1.15</td><td>✓</td></tr><tr><td>area average downward (CTL)</td><td>3.44</td><td>4.63</td><td>3.39</td><td>0.84</td><td>✓</td></tr><tr><td>area average upward (CC)</td><td>3.41 (+23.4%)</td><td>6.46 (+12.6%)</td><td>5.40 (+27.4%)</td><td>1.07 (-6.9%)</td><td>✓(✓)</td></tr><tr><td>area average downward (CC)</td><td>3.18 (-7.5%)</td><td>4.96 (+7.2%)</td><td>4.35 (+28.4%)</td><td>1.12 (+26.3%)</td><td>✓(-)</td></tr><tr><td>sensitivity (L<sup>2</sup>)</td><td>0.013</td><td>0.009</td><td>0.014</td><td>0.001</td><td>✓</td></tr><tr><td><strong>NO<sub>3</sub><sup>-</sup> concentration (mmol N m<sup>-3</sup>)</strong><br>(not shown)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>area average (CTL)</td><td>12.02</td><td>11.49</td><td>11.81</td><td>10.55</td><td>✓</td></tr><tr><td>area average (CC)</td><td>11.05 (-8.1%)</td><td>10.90 (-5.2%)</td><td>11.03 (-6.6%)</td><td>10.17 (-3.6%)</td><td>✓(✓)</td></tr><tr><td>sensitivity (L<sup>2</sup>)</td><td>1.00</td><td>0.65</td><td>0.86</td><td>0.64</td><td>✓</td></tr><tr><td><strong>NPP (mmol N m<sup>-2</sup> day<sup>-1</sup>)</strong><br>(Figs. 6 and 10)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>area average (CTL)</td><td>2.64</td><td>2.83</td><td>2.61</td><td>1.31</td><td>✓</td></tr><tr><td>area average (CC)</td><td>2.18 (-17.2%)</td><td>2.36 (-16.5%)</td><td>2.18 (-16.6%)</td><td>1.07 (-18.3%)</td><td>✓(-)</td></tr></tbody></table>

improves upon the use of a GM-based parameterisation as-is in the eddy-permitting regime (consistent with the results of Mak et al. 2023), there are cases where the sensitivities of the diagnostics are not necessarily improved by the use of the splitting algorithm, although there is also no strong evidence that the case with no parameterisation is better either. Nevertheless, the results are promising, supporting the conclusion that the use of a GM-based scheme in eddy-permitting is possible and desirable if the splitting procedure is employed. Practically, since the field splitting procedure via the application of the spatial filter is not performed every time-step (here it is performed every model day, on the assumption that the large-scale evolves on a slower time-scale), the computational costs are rather minimal, roughly around 5% additional cost for the present idealised model configuration at the eddy-permitting resolution (cf., Mak et al., 2023).

The present results do not invalidate the conclusions drawn in Ruan et al. (2023), which demonstrates that the use of the GM-based version of GEOMETRIC improves the modelled state and sensitivities as compared to a standard prescription of the GM-coefficient $\kappa_{gm}$ as a constant in a coarse resolution model. It is however largely true that the main factor leading to improvements in the modelled mean states appears to come from a model becoming eddy-permitting. Some part of the observed differences with coarse resolution models could be attributed to the presence of the non-trivial bathymetric slope or the absence of isoneutral diffusion (sample calculations not shown), but the results seem to suggest that the statement “a GM-based parameterisation turns an eddy-permitting model into an expensive coarse resolution model” that we hypothesised is overly simplified (e.g., the differences in the prescribed grid-scale viscosity). Using a parameterisation as-is appears

to degrade the representation of the mean states. For the present model, it seems to be possible to tune the parameters accordingly to reproduce a reasonable mean state (e.g., the GEOM calculation but with $\alpha = 0.025$), but the resulting model has an explicit eddy kinetic energy that is too low (the parameterisation here largely impacts the Western Boundary Current region; not shown), implying the variability has been affected. The splitting approach appears to be able to retain both of the desirable features of the explicit eddies and some action of the GM-based schemes in the present work and previous work of Mak et al. (2023), and displays aspects of scale-awareness that allow the parameterisation to be used across multiple grid resolutions without retuning. It is certainly true that the statistics of the resolved eddies could be different in eddy-permitting regimes relative to eddy-rich/resolving regimes, and quantifying the difference (as well as providing proposals for any fixes) should be considered in a future work.

The present work only investigates the use of splitting with a GM-based scheme, and does not consider more advanced procedures of tapering of $\kappa_{gm}$ as the shallow ocean is approached (e.g., Wei et al., 2022, 2024; Nummelin and Isachsen, 2024), although sample experiments with and without any tapering seem to make no significant quantitative difference to our results (not shown). The present work also does not employ isoneutral mixing (e.g., Redi, 1982; Griffies, 1998), which is known to modify the tracer transport rates and/or modify the ocean state (e.g., Jones and Abernathey, 2019; Holmes et al., 2022; Chouksey et al., 2022). One subtlety highlighted in Mak et al. (2023) is that extra mixing along the large-scale stratification profile may lead to significant diapycnal fluxes across the actual resolved stratification, and it is not clear whether the splitting approach
<caption>**Table 2**</caption>

Summary of diagnostics and their sensitivities for the set of calculations, all analysed within the analysis domain between $y = 25$ and $45^{\circ}$ N ($S_1$ and $S_3$ in the text) and $z = -700$ m, except for mixed layer depth, which is analysed north of $y = 45^{\circ}$ N in the deep water formation region. The bracketed numbers denote the percentage differences of the relevant diagnostic between the climate change and control scenario. $L^2$ sensitivity denotes the area-weighted average root-mean-square difference between the climate change and control scenario. A dash is given if there is no obvious evidence that SPLIT leads to an improvement over R4. Note that the values of the northward heat transport here are larger than those reported in Ruan et al. (2023), which was computed from an integral over the whole model depth, and was an inconsistent choice on our part.

<table><thead><tr><th>Diagnostic</th><th>R12 values</th><th>R4 values</th><th>SPLIT values</th><th>GEOM values</th><th>Improve by SPLIT</th></tr></thead><tbody><tr><td><strong>overturning circulation (Sv)</strong><br>(Figs. 4a-d and 9a-d)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>L<sup>2</sup> mismatch rel. R12 (CTL)</td><td>-</td><td>1.21</td><td>1.09</td><td>1.98</td><td>✓</td></tr><tr><td>L<sup>2</sup> mismatch rel. R12 (CC)</td><td>-</td><td>0.70</td><td>0.69</td><td>1.58</td><td>-</td></tr><tr><td><strong>northern mixed layer depth (m)</strong><br>(Figs. 4e-h and 9e-h)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>median (CTL)</td><td>504</td><td>765</td><td>664</td><td>311</td><td>✓</td></tr><tr><td>median (CC)</td><td>347 (-31.2%)</td><td>478 (-37.5%)</td><td>422 (-36.5%)</td><td>225 (-27.6%)</td><td>✓(✓)</td></tr><tr><td>quartile range (CTL)</td><td>232</td><td>439</td><td>401</td><td>197</td><td>✓</td></tr><tr><td>quartile range (CC)</td><td>111 (-52.4%)</td><td>199 (-54.6%)</td><td>172 (-57.1%)</td><td>128 (-34.8%)</td><td>✓(x)</td></tr><tr><td><strong>northward heat transport (PW)</strong><br>(Fig. 5a)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>area average (CTL)</td><td>3.89</td><td>7.29</td><td>4.41</td><td>-0.41</td><td>✓</td></tr><tr><td>area average (CC)</td><td>4.67 (+20.0%)</td><td>5.05 (-30.7%)</td><td>3.17 (-28.2%)</td><td>0.56 (-235.8%)</td><td>✓(-)</td></tr><tr><td>sensitivity (L<sup>2</sup>)</td><td>2.47</td><td>3.03</td><td>2.10</td><td>1.41</td><td>✓</td></tr><tr><td><strong>vertical heat transport (10<sup>-3</sup> PW)</strong><br>(Fig. 5b)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>area average upward (CTL)</td><td>2.76</td><td>5.74</td><td>4.24</td><td>1.15</td><td>✓</td></tr><tr><td>area average downward (CTL)</td><td>3.44</td><td>4.63</td><td>3.39</td><td>0.84</td><td>✓</td></tr><tr><td>area average upward (CC)</td><td>3.41 (+23.4%)</td><td>6.46 (+12.6%)</td><td>5.40 (+27.4%)</td><td>1.07 (-6.9%)</td><td>✓(✓)</td></tr><tr><td>area average downward (CC)</td><td>3.18 (-7.5%)</td><td>4.96 (+7.2%)</td><td>4.35 (+28.4%)</td><td>1.12 (+26.3%)</td><td>✓(-)</td></tr><tr><td>sensitivity (L<sup>2</sup>)</td><td>0.013</td><td>0.009</td><td>0.014</td><td>0.001</td><td>✓</td></tr><tr><td><strong>NO<sub>3</sub><sup>-</sup> concentration (mmol N m<sup>-3</sup>)</strong><br>(not shown)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>area average (CTL)</td><td>12.02</td><td>11.49</td><td>11.81</td><td>10.55</td><td>✓</td></tr><tr><td>area average (CC)</td><td>11.05 (-8.1%)</td><td>10.90 (-5.2%)</td><td>11.03 (-6.6%)</td><td>10.17 (-3.6%)</td><td>✓(✓)</td></tr><tr><td>sensitivity (L<sup>2</sup>)</td><td>1.00</td><td>0.65</td><td>0.86</td><td>0.64</td><td>✓</td></tr><tr><td><strong>NPP (mmol N m<sup>-2</sup> day<sup>-1</sup>)</strong><br>(Figs. 6 and 10)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>area average (CTL)</td><td>2.64</td><td>2.83</td><td>2.61</td><td>1.31</td><td>✓</td></tr><tr><td>area average (CC)</td><td>2.18 (-17.2%)</td><td>2.36 (-16.5%)</td><td>2.18 (-16.6%)</td><td>1.07 (-18.3%)</td><td>✓(-)</td></tr></tbody></table>

improves upon the use of a GM-based parameterisation as-is in the eddy-permitting regime (consistent with the results of Mak et al. 2023), there are cases where the sensitivities of the diagnostics are not necessarily improved by the use of the splitting algorithm, although there is also no strong evidence that the case with no parameterisation is better either. Nevertheless, the results are promising, supporting the conclusion that the use of a GM-based scheme in eddy-permitting is possible and desirable if the splitting procedure is employed. Practically, since the field splitting procedure via the application of the spatial filter is not performed every time-step (here it is performed every model day, on the assumption that the large-scale evolves on a slower time-scale), the computational costs are rather minimal, roughly around 5% additional cost for the present idealised model configuration at the eddy-permitting resolution (cf., Mak et al., 2023).

The present results do not invalidate the conclusions drawn in Ruan et al. (2023), which demonstrates that the use of the GM-based version of GEOMETRIC improves the modelled state and sensitivities as compared to a standard prescription of the GM-coefficient $\kappa_{gm}$ as a constant in a coarse resolution model. It is however largely true that the main factor leading to improvements in the modelled mean states appears to come from a model becoming eddy-permitting. Some part of the observed differences with coarse resolution models could be attributed to the presence of the non-trivial bathymetric slope or the absence of isoneutral diffusion (sample calculations not shown), but the results seem to suggest that the statement “a GM-based parameterisation turns an eddy-permitting model into an expensive coarse resolution model” that we hypothesised is overly simplified (e.g., the differences in the prescribed grid-scale viscosity). Using a parameterisation as-is appears

to degrade the representation of the mean states. For the present model, it seems to be possible to tune the parameters accordingly to reproduce a reasonable mean state (e.g., the GEOM calculation but with $\alpha = 0.025$), but the resulting model has an explicit eddy kinetic energy that is too low (the parameterisation here largely impacts the Western Boundary Current region; not shown), implying the variability has been affected. The splitting approach appears to be able to retain both of the desirable features of the explicit eddies and some action of the GM-based schemes in the present work and previous work of Mak et al. (2023), and displays aspects of scale-awareness that allow the parameterisation to be used across multiple grid resolutions without retuning. It is certainly true that the statistics of the resolved eddies could be different in eddy-permitting regimes relative to eddy-rich/resolving regimes, and quantifying the difference (as well as providing proposals for any fixes) should be considered in a future work.

The present work only investigates the use of splitting with a GM-based scheme, and does not consider more advanced procedures of tapering of $\kappa_{gm}$ as the shallow ocean is approached (e.g., Wei et al., 2022, 2024; Nummelin and Isachsen, 2024), although sample experiments with and without any tapering seem to make no significant quantitative difference to our results (not shown). The present work also does not employ isoneutral mixing (e.g., Redi, 1982; Griffies, 1998), which is known to modify the tracer transport rates and/or modify the ocean state (e.g., Jones and Abernathey, 2019; Holmes et al., 2022; Chouksey et al., 2022). One subtlety highlighted in Mak et al. (2023) is that extra mixing along the large-scale stratification profile may lead to significant diapycnal fluxes across the actual resolved stratification, and it is not clear whether the splitting approach
should be used in that context. The proposed safer option is that diffusion should remain along the *full* modelled isopycnal profile, either by saving extra variables during model run-time, or recomputing the isopycnal slopes. Sample diagnostics from calculations in the present model configuration at eddy-permitting resolution but with extra diffusion along the full and large-scale isopycnal profile are provided in *Appendix*. The use of isoneutral diffusion appear to lead to a damping of the explicit eddy activity, more so when the state-aware GEOMETRIC is employed. There appears to be a positive feedback loop, where isoneutral diffusion modifies the underlying tracer distribution and thus stratification, leading to changes in $\kappa_{gm}$ through GEOMETRIC, and modifying the stratification via the resulting eddy-induced velocity. The results presented in *Appendix* further illustrate the complexities of tuning and utilising parameterisations, and question a perceived view that increasing isoneutral diffusion increases the rate of ventilation but not necessarily the pathways, which is not true if state-aware GM-based schemes are utilised. A comprehensive exploration associated with the extra degrees of freedom and possible positive feedback loops from including isoneutral mixing is beyond the scope of the present paper, but is under investigation and will be reported in a future publication.

We have not considered in this work the inclusion of backscatter (e.g., Bachman, 2019; Jansen et al., 2019; Yankovsky et al., 2023), which could energise the Western Boundary Current, strengthen the too weak overturning in the subtropical gyre, and/or modify the tracer transport rates to further improve on both the physical and biogeochemical diagnostics of interest. The use of the splitting algorithm is not mutually exclusive of backscatter, and we refer the reader to the various subtleties one has to be aware of that has already been discussed in the work of Mak et al. (2023).

Probing and constraining the uncertainties of the splitting approach with GM-based GEOMETRIC scheme and its impacts on physical and biogeochemical responses is necessary with the increasing prevalence of eddy-permitting models. The assessments are essential for climate projections in realistic global configurations, initial conditions, atmospheric forcing, modelled biogeochemical processes, and so forth (e.g., Berthet et al., 2019; Swearer et al., 2019; Séférían et al., 2019; Couespel et al., 2024). The splitting approach provides one part of a solution to the problem of representing mesoscale eddy effects in numerical ocean modelling in eddy-permitting regime, and further investigations into other mesoscale eddy parameterisations are still necessary. An ongoing line of investigation relates to the use of the splitting procedure with a nonlinear equation of state and with backscatter in a global eddy-permitting model, and results from the associated research will be reported elsewhere in due course.

More generally, this work shows that our procedure of interest leads to an improvement in the modelled biogeochemical response in a physically consistent and expected way. We advocate that similar assessments for both the modelled mean state and its sensitivities (to forcing scenarios, free parameters, or otherwise) being performed in relation to other parameterisation approaches, be they deterministic, stochastic and/or data-driven. As highlighted in Ruan et al. (2023), a reasonable biogeochemical response could arise from a physically inconsistent response. Ultimately one should be aiming at procedures that get the ‘right’ answer for the ‘right’ reasons, and such assessments provide evidence in support of a procedure’s soundness and/or robustness in other regimes.

### CRediT authorship contribution statement

**X. Ruan:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation. **D. Couespel:** Writing – review & editing, Writing – original draft, Software, Methodology, Conceptualization. **M. Lévy:** Writing – review & editing, Writing – original draft, Methodology. **J. Li:** Writing – review & editing, Writing – original draft, Supervision. **J. Mak:** Writing – review & editing, Writing – original draft,

Visualization, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization. **Y. Wang:** Writing – review & editing, Writing – original draft, Supervision.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

Data is shared and published at http://dx.doi.org/10.5281/zenodo.11498192.

### Acknowledgements

This research was funded by both the RGC General Research Fund 16304021 and the Center for Ocean Research in Hong Kong and Macau, a joint research centre for ocean research between Laoshan Laboratory and Hong Kong University of Science and Technology. The numerical model modifications, analyses code and sample model data may be found on Zenodo at http://dx.doi.org/10.5281/zenodo.11498192. We thank the two anonymous referees for their supporting comments that have lead to an improvement in the presentation of the article.

### Appendix. Inclusion of isoneutral diffusion

In this section, we provide sample numerical results for the R4 and SPLIT calculations employing isoneutral diffusion; analogous results for GEOM have been omitted, since the results are largely similar to that reported in the Appendix of Ruan et al. (2023). When a harmonic tracer diffusion along a isoneutral direction is switched on, the biharmonic horizontal tracer diffusion is switched off; the harmonic operator is expected to largely supersede the action of the biharmonic operator since it damps over a broader length-scale (cf. in Fourier space, these would correspond to a damping of $-k^2$ and $-k^4$ where $k$ is a representative wavenumber). A diffusivity of $\kappa_{iso} = 500 \text{ m}^2 \text{ s}^{-1}$ is utilised which, while possibly on the slightly large side, demonstrates on why the use of isoneutral diffusion together with a state-aware GM-based parameterisation requires more care than is perhaps acknowledged in the literature.

As noted in Mak et al. (2023), diffusing along the isoneutral direction associated with the large-scale isopycnals might lead to uncontrolled dianeutral fluxes, and a safer option is to diffuse along the isoneutral direction associated with the *full* isopycnals. For completeness, however, we tried both approaches, with experiments termed R4(Redi_full) and SPLIT(Redi_full), which employs isoneutral diffusion along the full isopycnals, and SPLIT(Redi_large), which has isoneutral diffusion along the large-scale isopycnals computed from the splitting algorithm. To get a sense of the immediate impacts of employing isoneutral diffusion, we show in Fig. A.12a, b, d, e, g the surface relative vorticity for the relevant calculations without and with isoneutral diffusion. The presence of isoneutral diffusion in the R4 cases (panels a and d) leads to some damping of the explicit fluctuations (although possibly not to the same extent as that reported in Mak et al. 2023; see their Fig. 2b). The damping might be expected, since isoneutral diffusion would erode tracer gradients, which by geostrophic balance would have an impact on the flow field of the baroclinic eddies, such as that seen in the relative vorticity field. A more dramatic damping is seen when GEOM is utilised, regardless of whether isoneutral diffusion is along the full isopycnals (panel e) or the large-scale isopycnals (panel g), relative to the case with no isopycnal diffusion (panel b). SPLIT(Redi_full) appears to experience the largest damping, some of which can presumably be attributed to the eddy-induced velocity via increases in the resulting
should be used in that context. The proposed safer option is that diffusion should remain along the *full* modelled isopycnal profile, either by saving extra variables during model run-time, or recomputing the isopycnal slopes. Sample diagnostics from calculations in the present model configuration at eddy-permitting resolution but with extra diffusion along the full and large-scale isopycnal profile are provided in *Appendix*. The use of isoneutral diffusion appear to lead to a damping of the explicit eddy activity, more so when the state-aware GEOMETRIC is employed. There appears to be a positive feedback loop, where isoneutral diffusion modifies the underlying tracer distribution and thus stratification, leading to changes in $\kappa_{gm}$ through GEOMETRIC, and modifying the stratification via the resulting eddy-induced velocity. The results presented in *Appendix* further illustrate the complexities of tuning and utilising parameterisations, and question a perceived view that increasing isoneutral diffusion increases the rate of ventilation but not necessarily the pathways, which is not true if state-aware GM-based schemes are utilised. A comprehensive exploration associated with the extra degrees of freedom and possible positive feedback loops from including isoneutral mixing is beyond the scope of the present paper, but is under investigation and will be reported in a future publication.

We have not considered in this work the inclusion of backscatter (e.g., Bachman, 2019; Jansen et al., 2019; Yankovsky et al., 2023), which could energise the Western Boundary Current, strengthen the too weak overturning in the subtropical gyre, and/or modify the tracer transport rates to further improve on both the physical and biogeochemical diagnostics of interest. The use of the splitting algorithm is not mutually exclusive of backscatter, and we refer the reader to the various subtleties one has to be aware of that has already been discussed in the work of Mak et al. (2023).

Probing and constraining the uncertainties of the splitting approach with GM-based GEOMETRIC scheme and its impacts on physical and biogeochemical responses is necessary with the increasing prevalence of eddy-permitting models. The assessments are essential for climate projections in realistic global configurations, initial conditions, atmospheric forcing, modelled biogeochemical processes, and so forth (e.g., Berthet et al., 2019; Swearer et al., 2019; Séférían et al., 2019; Couespel et al., 2024). The splitting approach provides one part of a solution to the problem of representing mesoscale eddy effects in numerical ocean modelling in eddy-permitting regime, and further investigations into other mesoscale eddy parameterisations are still necessary. An ongoing line of investigation relates to the use of the splitting procedure with a nonlinear equation of state and with backscatter in a global eddy-permitting model, and results from the associated research will be reported elsewhere in due course.

More generally, this work shows that our procedure of interest leads to an improvement in the modelled biogeochemical response in a physically consistent and expected way. We advocate that similar assessments for both the modelled mean state and its sensitivities (to forcing scenarios, free parameters, or otherwise) being performed in relation to other parameterisation approaches, be they deterministic, stochastic and/or data-driven. As highlighted in Ruan et al. (2023), a reasonable biogeochemical response could arise from a physically inconsistent response. Ultimately one should be aiming at procedures that get the ‘right’ answer for the ‘right’ reasons, and such assessments provide evidence in support of a procedure’s soundness and/or robustness in other regimes.

### CRediT authorship contribution statement

**X. Ruan:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation. **D. Couespel:** Writing – review & editing, Writing – original draft, Software, Methodology, Conceptualization. **M. Lévy:** Writing – review & editing, Writing – original draft, Methodology. **J. Li:** Writing – review & editing, Writing – original draft, Supervision. **J. Mak:** Writing – review & editing, Writing – original draft,

Visualization, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization. **Y. Wang:** Writing – review & editing, Writing – original draft, Supervision.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

Data is shared and published at http://dx.doi.org/10.5281/zenodo.11498192.

### Acknowledgements

This research was funded by both the RGC General Research Fund 16304021 and the Center for Ocean Research in Hong Kong and Macau, a joint research centre for ocean research between Laoshan Laboratory and Hong Kong University of Science and Technology. The numerical model modifications, analyses code and sample model data may be found on Zenodo at http://dx.doi.org/10.5281/zenodo.11498192. We thank the two anonymous referees for their supporting comments that have lead to an improvement in the presentation of the article.

## Appendix. Inclusion of isoneutral diffusion

In this section, we provide sample numerical results for the R4 and SPLIT calculations employing isoneutral diffusion; analogous results for GEOM have been omitted, since the results are largely similar to that reported in the Appendix of Ruan et al. (2023). When a harmonic tracer diffusion along a isoneutral direction is switched on, the biharmonic horizontal tracer diffusion is switched off; the harmonic operator is expected to largely supersede the action of the biharmonic operator since it damps over a broader length-scale (cf. in Fourier space, these would correspond to a damping of $-k^2$ and $-k^4$ where $k$ is a representative wavenumber). A diffusivity of $\kappa_{iso} = 500 \text{ m}^2 \text{ s}^{-1}$ is utilised which, while possibly on the slightly large side, demonstrates on why the use of isoneutral diffusion together with a state-aware GM-based parameterisation requires more care than is perhaps acknowledged in the literature.

As noted in Mak et al. (2023), diffusing along the isoneutral direction associated with the large-scale isopycnals might lead to uncontrolled dianeutral fluxes, and a safer option is to diffuse along the isoneutral direction associated with the *full* isopycnals. For completeness, however, we tried both approaches, with experiments termed R4(Redi_full) and SPLIT(Redi_full), which employs isoneutral diffusion along the full isopycnals, and SPLIT(Redi_large), which has isoneutral diffusion along the large-scale isopycnals computed from the splitting algorithm. To get a sense of the immediate impacts of employing isoneutral diffusion, we show in Fig. A.12a, b, d, e, g the surface relative vorticity for the relevant calculations without and with isoneutral diffusion. The presence of isoneutral diffusion in the R4 cases (panels a and d) leads to some damping of the explicit fluctuations (although possibly not to the same extent as that reported in Mak et al. 2023; see their Fig. 2b). The damping might be expected, since isoneutral diffusion would erode tracer gradients, which by geostrophic balance would have an impact on the flow field of the baroclinic eddies, such as that seen in the relative vorticity field. A more dramatic damping is seen when GEOM is utilised, regardless of whether isoneutral diffusion is along the full isopycnals (panel e) or the large-scale isopycnals (panel g), relative to the case with no isopycnal diffusion (panel b). SPLIT(Redi_full) appears to experience the largest damping, some of which can presumably be attributed to the eddy-induced velocity via increases in the resulting
<figure><img src="image_13.png" /><figcaption>**Fig. A.12.** Snapshots of the surface relative vorticity (in units of the planetary vorticity $f_0$) for (a, d) the R4 cases, and (b, e, g) the SPLIT cases, without isoneutral diffusion, with isoneutral diffusion along the full isopycnal slopes, and with isoneutral along the large-scale isopycnal slopes. We show in (c, f, h) the corresponding $\kappa_{gm}(x, y, t)$ values (in units of $m^2 s^{-1}$) for the SPLIT cases. Where isoneutral diffusion is active, the value of the diffusivity is $\kappa_{gm} = 500 \ m^2 s^{-1}$.</figcaption></figure>

$\kappa_{gm}$ (panel f). There is still noticeable damping in SPLIT(Redi_large) (panel g) even if the changes in $\kappa_{gm}$ are rather mild (panel h), providing extra evidence that the isoneutral diffusion has an effect on the explicit variability. The associated damping can also be seen in the time-averaged sense by measures of the domain-integrated explicit eddy kinetic energy for example (not shown): the decrease in domain-integrated eddy kinetic energy is found to be of a larger percentage in the corresponding SPLIT experiments compared to the corresponding R4 experiments.

This change in $\kappa_{gm}$ we argue to result from isoneutral diffusion modifying the tracer distribution and the stratification, leading to changes in $\kappa_{gm}$ via the GEOMETRIC parameterisation (through both the prescription of $\kappa_{gm}$ in Eq. (5) and via the eddy energy budget in Eq. (6)), which further leads to changes in the stratification and modifying $\kappa_{gm}$. Such a positive feedback loop is absent in cases where there is no GM-based parameterisation employed (as in R4 here), or in cases where $\kappa_{gm}$ is fixed (as seen in Appendix of Ruan et al. 2023); the splitting approach reduces the degree of (but does not remove) this positive feedback loop.

Table A.3 shows the relevant metrics diagnosed from the R4 and SPLIT calculations with isoneutral diffusion active, to be compared with the corresponding values in Table 2. As a summary, it may be seen

that the R4 calculation with isoneutral diffusion remain very similar to that without, with relatively minor *increases* in the relevant transports (meridional and vertical heat transports, increased nitrate concentration), resulting in marginally larger NPP values. The same cannot be said of the SPLIT calculation with isoneutral diffusion. SPLIT(Redi_full) experiences large shoaling of the mixed layer depths (the median decreases by about 200 meters in the control case, and about 100 meters in the climate change case), and a significant decrease in both the meridional and vertical heat transports, which are symptoms of a substantially reduced meridional overturning circulation consistent via changes in the stratification from increases in the value of $\kappa_{gm}$ (cf. observations in Ruan et al. 2023). Curiously, the resulting nitrate concentration increases somewhat. However, the corresponding NPP value noticeably decreases, and the decrease in NPP with increasing $\kappa_{gm}$ is consistent with the results from Ruan et al. (2023). On the other hand, the metrics associated with SPLIT(Redi_large) are not entirely unlike that of R4(Redi_full), with an overly deep mixed layer, rather large transport, and a rather large NPP. This could have arisen from the additional transport due to the isoneutral diffusion, but also possibly from the expected spurious diapycnal mixing arising from the isoneutral diffusion but along the direction of the large-scale rather than full
<figure><img src="image_14.png" /><figcaption>**Fig. A.12.** Snapshots of the surface relative vorticity (in units of the planetary vorticity $f_0$) for (a, d) the R4 cases, and (b, e, g) the SPLIT cases, without isoneutral diffusion, with isoneutral diffusion along the full isopycnal slopes, and with isoneutral along the large-scale isopycnal slopes. We show in (c, f, h) the corresponding $\kappa_{gm}(x, y, t)$ values (in units of $m^2 s^{-1}$) for the SPLIT cases. Where isoneutral diffusion is active, the value of the diffusivity is $\kappa_{gm} = 500 \ m^2 s^{-1}$.</figcaption></figure>

$\kappa_{gm}$ (panel f). There is still noticeable damping in SPLIT(Redi_large) (panel g) even if the changes in $\kappa_{gm}$ are rather mild (panel h), providing extra evidence that the isoneutral diffusion has an effect on the explicit variability. The associated damping can also be seen in the time-averaged sense by measures of the domain-integrated explicit eddy kinetic energy for example (not shown): the decrease in domain-integrated eddy kinetic energy is found to be of a larger percentage in the corresponding SPLIT experiments compared to the corresponding R4 experiments.

This change in $\kappa_{gm}$ we argue to result from isoneutral diffusion modifying the tracer distribution and the stratification, leading to changes in $\kappa_{gm}$ via the GEOMETRIC parameterisation (through both the prescription of $\kappa_{gm}$ in Eq. (5) and via the eddy energy budget in Eq. (6)), which further leads to changes in the stratification and modifying $\kappa_{gm}$. Such a positive feedback loop is absent in cases where there is no GM-based parameterisation employed (as in R4 here), or in cases where $\kappa_{gm}$ is fixed (as seen in Appendix of Ruan et al. 2023); the splitting approach reduces the degree of (but does not remove) this positive feedback loop.

Table A.3 shows the relevant metrics diagnosed from the R4 and SPLIT calculations with isoneutral diffusion active, to be compared with the corresponding values in Table 2. As a summary, it may be seen

that the R4 calculation with isoneutral diffusion remain very similar to that without, with relatively minor *increases* in the relevant transports (meridional and vertical heat transports, increased nitrate concentration), resulting in marginally larger NPP values. The same cannot be said of the SPLIT calculation with isoneutral diffusion. SPLIT(Redi_full) experiences large shoaling of the mixed layer depths (the median decreases by about 200 meters in the control case, and about 100 meters in the climate change case), and a significant decrease in both the meridional and vertical heat transports, which are symptoms of a substantially reduced meridional overturning circulation consistent via changes in the stratification from increases in the value of $\kappa_{gm}$ (cf. observations in Ruan et al. 2023). Curiously, the resulting nitrate concentration increases somewhat. However, the corresponding NPP value noticeably decreases, and the decrease in NPP with increasing $\kappa_{gm}$ is consistent with the results from Ruan et al. (2023). On the other hand, the metrics associated with SPLIT(Redi_large) are not entirely unlike that of R4(Redi_full), with an overly deep mixed layer, rather large transport, and a rather large NPP. This could have arisen from the additional transport due to the isoneutral diffusion, but also possibly from the expected spurious diapycnal mixing arising from the isoneutral diffusion but along the direction of the large-scale rather than full
<caption>Table A.3</caption>

Summary of diagnostics and their sensitivities for calculation with isoneutral diffusion active in R4, and active on the full-scale in SPLIT, and active on the large-scale in the SPLIT with codes modified ($k_{iso} = 500 \text{ m}^2 \text{ s}^{-1}$ in all of these calculations), to be compared to values in Table 2. All metrics were diagnosed within the analysis domain between $y = 25$ and $45^\circ \text{ N}$ ($S_1$ and $S_3$ in the text) and $z = -700 \text{ m}$, except for mixed layer depth, which is analysed north of $y = 45^\circ \text{ N}$ in the deep water formation region. The bracketed numbers denote the percentage differences of the relevant diagnostic between the climate change and control scenario. $L^2$ sensitivity denotes the area-weighted average root-mean-square difference between the climate change and control scenario.

<table><thead><tr><th>Diagnostic</th><th>R4(Redi_full) values</th><th>SPLIT(Redi_full) values</th><th>SPLIT(Redi_large) values</th></tr></thead><tbody><tr><td><strong>overturning circulation (Sv)</strong></td><td></td><td></td><td></td></tr><tr><td>L<sup>2</sup> mismatch rel. R12 (CTL)</td><td>1.25</td><td>1.36</td><td>1.71</td></tr><tr><td>L<sup>2</sup> mismatch rel. R12 (CC)</td><td>0.65</td><td>1.08</td><td>0.98</td></tr><tr><td><strong>northern mixed layer depth (m)</strong></td><td></td><td></td><td></td></tr><tr><td>median (CTL)</td><td>773</td><td>479</td><td>930</td></tr><tr><td>median (CC)</td><td>461 (-40.4%)</td><td>292 (-39.0%)</td><td>494 (-46.9%)</td></tr><tr><td>quartile range (CTL)</td><td>492</td><td>332</td><td>551</td></tr><tr><td>quartile range (CC)</td><td>211 (-57.0%)</td><td>174 (-47.6%)</td><td>239 (-56.6%)</td></tr><tr><td><strong>northward heat transport (PW)</strong></td><td></td><td></td><td></td></tr><tr><td>area average (CTL)</td><td>7.49</td><td>1.08</td><td>6.57</td></tr><tr><td>area average (CC)</td><td>5.20 (-30.6%)</td><td>1.63 (+50.9%)</td><td>3.54 (-46.1%)</td></tr><tr><td>sensitivity (L<sup>2</sup>)</td><td>2.97</td><td>1.66</td><td>3.57</td></tr><tr><td><strong>vertical heat transport (10<sup>-3</sup> PW)</strong></td><td></td><td></td><td></td></tr><tr><td>area average upward (CTL)</td><td>6.03</td><td>1.86</td><td>4.36</td></tr><tr><td>area average downward (CTL)</td><td>4.80</td><td>1.54</td><td>2.62</td></tr><tr><td>area average upward (CC)</td><td>6.17 (+2.3%)</td><td>1.49 (-20.1%)</td><td>5.59 (+28.2%)</td></tr><tr><td>area average downward (CC)</td><td>5.05 (+5.2%)</td><td>1.38 (-10.4%)</td><td>4.28 (+63.1%)</td></tr><tr><td>sensitivity (L<sup>2</sup>)</td><td>0.006</td><td>0.003</td><td>0.006</td></tr><tr><td><strong>NO<sub>3</sub><sup>-</sup> concentration (mmol N m<sup>-3</sup>)</strong></td><td></td><td></td><td></td></tr><tr><td>area average (CTL)</td><td>11.61</td><td>12.56</td><td>11.82</td></tr><tr><td>area average (CC)</td><td>11.01 (-5.1%)</td><td>11.76 (-6.4%)</td><td>11.30 (-4.4%)</td></tr><tr><td>sensitivity (L<sup>2</sup>)</td><td>0.65</td><td>1.03</td><td>0.58</td></tr><tr><td><strong>NPP (mmol N m<sup>-2</sup> day<sup>-1</sup>)</strong></td><td></td><td></td><td></td></tr><tr><td>area average (CTL)</td><td>2.88</td><td>2.13</td><td>2.89</td></tr><tr><td>area average (CC)</td><td>2.42 (-15.9%)</td><td>1.62 (-23.9%)</td><td>2.50 (-13.6%)</td></tr></tbody></table>

isopycnals. The latter is somewhat harder to quantify and is beyond the scope of the present work, although there are frameworks for doing so (e.g., Lee et al., 2002; Megann, 2018).

The results here highlight further complexities in utilising isoneutral diffusion and eddy-induced advection parameterisations together, particularly when state-aware parameterisations are used, because of possible feedback loops with the increased complexity of the parameterisations. In the present model, inclusion of isoneutral diffusion has rather weak effect on the R4 calculations, and there appears to be a positive feedback loop present when a state-aware parameterisation for the eddy-induced advection such as GEOMETRIC is utilised. Further work is still needed on methodologies for tuning strategies for the parameterisations, or possibly on a parameterisation that unifies the two separate but dynamically related processes, which are really manifestations of the same underlying baroclinic turbulence.

## References

Bachman, S.D., 2019. The GM+E closure: A framework for coupling backscatter with the Gent and McWilliams parameterization. Ocean Model. 136, 85-106.

Beech, N., Rackow, T., Semmler, T., Danilov, S., Wang, Q., Jung, T., 2022. Long-term evolution of ocean eddy activity in a warming world. Nat. Clim. Change 12 (10), 910-917.

Berthet, S., Séférían, R., Bricaud, C., Chevallier, M., Voldoire, A., Ethé, C., 2019. Evaluation of an online grid-coarsening algorithm in a global eddy-admitting ocean biogeochemical model. J. Adv. Model. Earth Syst. 11, 1-25.

Chang, C., Adcroft, A., Zanna, L., Hallberg, R., Griffies, S.M., 2023. Remote versus local impacts of energy backscatter on the North Atlantic SST biases in a global ocean model. Geophys. Res. Lett. 50, e2023GL105757.

Charney, J.G., 1971. Geostrophic turbulence. J. Atmos. Sci. 28, 1087-1095.

Chelton, D.B., Schlax, M.G., Samelson, R.M., 2011. Global observations of nonlinear mesoscale eddies. Prog. Oceanog. 91, 167-216.

Chouksey, A., Grisel, A., Chouksey, M., Eden, C., 2022. Changes in global ocean circulation due to isopycnal diffusion. J. Phys. Oceanogr. 52, 2219-2235.

Couespel, D., Lévy, M., Bopp, L., 2021. Oceanic primary production decline halved in eddy-resolving simulations of global warming. Biogeosciences 18, 4321-4349.

Couespel, D., Lévy, M., Bopp, L., 2024. Stronger oceanic CO<sub>2</sub> sink in eddy-resolving simulations of global warming. Geophys. Res. Lett. 51, e2023GL106172.

Eden, C., Greatbatch, R.J., 2008. Towards a mesoscale eddy closure. Ocean Model. 20, 223-239.

Fox-Kemper, B., Adcroft, A.J., Böning, C.W., Chassignet, E.P., Curchitser, E.N., Danabasoglu, G., Eden, C., England, M.H., Gerdes, R., Greatbatch, R.J., Griffies, S.M., Hallberg, R.W., Hanert, E., Heimbach, P., Hewitt, H.T., Hill, C.N., Komuro, Y., Legg, S., Le Sommer, J., Masina, S., Marsland, S.J., Penny, S.G., Qiao, F., Ringler, T.D., Treguier, A.M., Tsujino, H., Uotila, P., Yeager, S.G., 2019. Challenges and prospects in ocean circulation models. Front. Mar. Sci 6, 65.

Gaspar, P., Grégoris, Y., Lefevre, J., 1990. A simple eddy kinetic energy model for simulations of the oceanic vertical mixing: Tests at station papa and long-term upper ocean study site. J. Geophys. Res. 95, 16179-16193.

Gent, P.R., McWilliams, J.C., 1990. Isopycnal mixing in ocean circulation models. J. Phys. Oceanogr. 20, 150-155.

Gent, P.R., Willebrand, J., McDougall, T.J., McWilliams, J.C., 1995. Parameterizing eddy-induced tracer transports in ocean circulation models. J. Phys. Oceanogr. 25, 463-474.

Griffies, S.M., 1998. The Gent-McWilliams skew flux. J. Phys. Oceanogr. 28, 831-841.

Griffies, S.M., Winton, M., Anderson, W.G., Benson, R., Delworth, T.L., Dufour, C.O., Dunne, J.P., Goddard, P., Morrison, A.K., Rosati, A., Wittenberg, A.T., Yin, J., Zhang, R., 2015. Impacts on ocean heat from transient mesoscale eddies in a hierarchy of climate models. J. Clim. 28 (3), 952-977.

Grooms, I., 2015. A computational study of turbulent kinetic energy transport in barotropic turbulence on the f-plane. Phys. Fluids 27, 101701.

Hallberg, R., 2013. Using a resolution function to regulate parameterizations of oceanic mesoscale eddy effects. Ocean Model. 72, 92-103.

Hewitt, H.T., Fox-Kemper, B., Pearson, B., Roberts, M., Klocke, D., 2022. The small scales of the ocean may hold the key to surprises. Nature Clim. Change 12, 496-499.

Hewitt, H.T., Roberts, M., Mathiot, P., Biastoch, A., Blockley, E., Chassignet, E.P., Fox-Kemper, B., Hyder, P., Marshall, D.P., Popova, E., Treguier, A., Zanna, L., Yool, A., Yu, Y., Beadling, R., Bell, M.J., Kuhlbrodt, T., Arsouze, T., Bellucci, A., Castruccio, F., Gan, B., Putrasahan, D., Roberts, C.D., Van Roekel, L., Zhang, Q., 2020. Resolving and parameterising the ocean mesoscale in earth system models. Curr. Clim. Change Rep. 6, 137-152.

Hoegh-Guldberg, O.R., Cai, R., Poloczanska, E.S., Brewer, P.G., Sundby, S., Hilmi, K., Fabry, V.J., Jung, S., 2014. The ocean. In: Barros, V.R., Field, C.B., Dokken, D.J., Mastrandrea, M.D., Mach, K.J., Bilir, T.E., Chatterjee, M., Ebi, K.L., Estrada, Y.O., Genova, R.C., Girma, B., Kissel, A.N., MacCracken, S., Mastrandrea, P.R., White, L.
<caption>Table A.3</caption>

Summary of diagnostics and their sensitivities for calculation with isoneutral diffusion active in R4, and active on the full-scale in SPLIT, and active on the large-scale in the SPLIT with codes modified ($k_{iso} = 500 \text{ m}^2 \text{ s}^{-1}$ in all of these calculations), to be compared to values in Table 2. All metrics were diagnosed within the analysis domain between $y = 25$ and $45^\circ \text{ N}$ ($S_1$ and $S_3$ in the text) and $z = -700 \text{ m}$, except for mixed layer depth, which is analysed north of $y = 45^\circ \text{ N}$ in the deep water formation region. The bracketed numbers denote the percentage differences of the relevant diagnostic between the climate change and control scenario. $L^2$ sensitivity denotes the area-weighted average root-mean-square difference between the climate change and control scenario.

<table><thead><tr><th>Diagnostic</th><th>R4(Redi_full) values</th><th>SPLIT(Redi_full) values</th><th>SPLIT(Redi_large) values</th></tr></thead><tbody><tr><td><strong>overturning circulation (Sv)</strong></td><td></td><td></td><td></td></tr><tr><td>L<sup>2</sup> mismatch rel. R12 (CTL)</td><td>1.25</td><td>1.36</td><td>1.71</td></tr><tr><td>L<sup>2</sup> mismatch rel. R12 (CC)</td><td>0.65</td><td>1.08</td><td>0.98</td></tr><tr><td><strong>northern mixed layer depth (m)</strong></td><td></td><td></td><td></td></tr><tr><td>median (CTL)</td><td>773</td><td>479</td><td>930</td></tr><tr><td>median (CC)</td><td>461 (-40.4%)</td><td>292 (-39.0%)</td><td>494 (-46.9%)</td></tr><tr><td>quartile range (CTL)</td><td>492</td><td>332</td><td>551</td></tr><tr><td>quartile range (CC)</td><td>211 (-57.0%)</td><td>174 (-47.6%)</td><td>239 (-56.6%)</td></tr><tr><td><strong>northward heat transport (PW)</strong></td><td></td><td></td><td></td></tr><tr><td>area average (CTL)</td><td>7.49</td><td>1.08</td><td>6.57</td></tr><tr><td>area average (CC)</td><td>5.20 (-30.6%)</td><td>1.63 (+50.9%)</td><td>3.54 (-46.1%)</td></tr><tr><td>sensitivity (L<sup>2</sup>)</td><td>2.97</td><td>1.66</td><td>3.57</td></tr><tr><td><strong>vertical heat transport (10<sup>-3</sup> PW)</strong></td><td></td><td></td><td></td></tr><tr><td>area average upward (CTL)</td><td>6.03</td><td>1.86</td><td>4.36</td></tr><tr><td>area average downward (CTL)</td><td>4.80</td><td>1.54</td><td>2.62</td></tr><tr><td>area average upward (CC)</td><td>6.17 (+2.3%)</td><td>1.49 (-20.1%)</td><td>5.59 (+28.2%)</td></tr><tr><td>area average downward (CC)</td><td>5.05 (+5.2%)</td><td>1.38 (-10.4%)</td><td>4.28 (+63.1%)</td></tr><tr><td>sensitivity (L<sup>2</sup>)</td><td>0.006</td><td>0.003</td><td>0.006</td></tr><tr><td><strong>NO<sub>3</sub><sup>-</sup> concentration (mmol N m<sup>-3</sup>)</strong></td><td></td><td></td><td></td></tr><tr><td>area average (CTL)</td><td>11.61</td><td>12.56</td><td>11.82</td></tr><tr><td>area average (CC)</td><td>11.01 (-5.1%)</td><td>11.76 (-6.4%)</td><td>11.30 (-4.4%)</td></tr><tr><td>sensitivity (L<sup>2</sup>)</td><td>0.65</td><td>1.03</td><td>0.58</td></tr><tr><td><strong>NPP (mmol N m<sup>-2</sup> day<sup>-1</sup>)</strong></td><td></td><td></td><td></td></tr><tr><td>area average (CTL)</td><td>2.88</td><td>2.13</td><td>2.89</td></tr><tr><td>area average (CC)</td><td>2.42 (-15.9%)</td><td>1.62 (-23.9%)</td><td>2.50 (-13.6%)</td></tr></tbody></table>

isopycnals. The latter is somewhat harder to quantify and is beyond the scope of the present work, although there are frameworks for doing so (e.g., Lee et al., 2002; Megann, 2018).

The results here highlight further complexities in utilising isoneutral diffusion and eddy-induced advection parameterisations together, particularly when state-aware parameterisations are used, because of possible feedback loops with the increased complexity of the parameterisations. In the present model, inclusion of isoneutral diffusion has rather weak effect on the R4 calculations, and there appears to be a positive feedback loop present when a state-aware parameterisation for the eddy-induced advection such as GEOMETRIC is utilised. Further work is still needed on methodologies for tuning strategies for the parameterisations, or possibly on a parameterisation that unifies the two separate but dynamically related processes, which are really manifestations of the same underlying baroclinic turbulence.

## References

Bachman, S.D., 2019. The GM+E closure: A framework for coupling backscatter with the Gent and McWilliams parameterization. Ocean Model. 136, 85-106.

Beech, N., Rackow, T., Semmler, T., Danilov, S., Wang, Q., Jung, T., 2022. Long-term evolution of ocean eddy activity in a warming world. Nat. Clim. Change 12 (10), 910-917.

Berthet, S., Séférían, R., Bricaud, C., Chevallier, M., Voldoire, A., Ethé, C., 2019. Evaluation of an online grid-coarsening algorithm in a global eddy-admitting ocean biogeochemical model. J. Adv. Model. Earth Syst. 11, 1-25.

Chang, C., Adcroft, A., Zanna, L., Hallberg, R., Griffies, S.M., 2023. Remote versus local impacts of energy backscatter on the North Atlantic SST biases in a global ocean model. Geophys. Res. Lett. 50, e2023GL105757.

Charney, J.G., 1971. Geostrophic turbulence. J. Atmos. Sci. 28, 1087-1095.

Chelton, D.B., Schlax, M.G., Samelson, R.M., 2011. Global observations of nonlinear mesoscale eddies. Prog. Oceanog. 91, 167-216.

Chouksey, A., Grisel, A., Chouksey, M., Eden, C., 2022. Changes in global ocean circulation due to isopycnal diffusion. J. Phys. Oceanogr. 52, 2219-2235.

Couespel, D., Lévy, M., Bopp, L., 2021. Oceanic primary production decline halved in eddy-resolving simulations of global warming. Biogeosciences 18, 4321-4349.

Couespel, D., Lévy, M., Bopp, L., 2024. Stronger oceanic CO<sub>2</sub> sink in eddy-resolving simulations of global warming. Geophys. Res. Lett. 51, e2023GL106172.

Eden, C., Greatbatch, R.J., 2008. Towards a mesoscale eddy closure. Ocean Model. 20, 223-239.

Fox-Kemper, B., Adcroft, A.J., Böning, C.W., Chassignet, E.P., Curchitser, E.N., Danabasoglu, G., Eden, C., England, M.H., Gerdes, R., Greatbatch, R.J., Griffies, S.M., Hallberg, R.W., Hanert, E., Heimbach, P., Hewitt, H.T., Hill, C.N., Komuro, Y., Legg, S., Le Sommer, J., Masina, S., Marsland, S.J., Penny, S.G., Qiao, F., Ringler, T.D., Treguier, A.M., Tsujino, H., Uotila, P., Yeager, S.G., 2019. Challenges and prospects in ocean circulation models. Front. Mar. Sci 6, 65.

Gaspar, P., Grégoris, Y., Lefevre, J., 1990. A simple eddy kinetic energy model for simulations of the oceanic vertical mixing: Tests at station papa and long-term upper ocean study site. J. Geophys. Res. 95, 16179-16193.

Gent, P.R., McWilliams, J.C., 1990. Isopycnal mixing in ocean circulation models. J. Phys. Oceanogr. 20, 150-155.

Gent, P.R., Willebrand, J., McDougall, T.J., McWilliams, J.C., 1995. Parameterizing eddy-induced tracer transports in ocean circulation models. J. Phys. Oceanogr. 25, 463-474.

Griffies, S.M., 1998. The Gent-McWilliams skew flux. J. Phys. Oceanogr. 28, 831-841.

Griffies, S.M., Winton, M., Anderson, W.G., Benson, R., Delworth, T.L., Dufour, C.O., Dunne, J.P., Goddard, P., Morrison, A.K., Rosati, A., Wittenberg, A.T., Yin, J., Zhang, R., 2015. Impacts on ocean heat from transient mesoscale eddies in a hierarchy of climate models. J. Clim. 28 (3), 952-977.

Grooms, I., 2015. A computational study of turbulent kinetic energy transport in barotropic turbulence on the f-plane. Phys. Fluids 27, 101701.

Hallberg, R., 2013. Using a resolution function to regulate parameterizations of oceanic mesoscale eddy effects. Ocean Model. 72, 92-103.

Hewitt, H.T., Fox-Kemper, B., Pearson, B., Roberts, M., Klocke, D., 2022. The small scales of the ocean may hold the key to surprises. Nature Clim. Change 12, 496-499.

Hewitt, H.T., Roberts, M., Mathiot, P., Biastoch, A., Blockley, E., Chassignet, E.P., Fox-Kemper, B., Hyder, P., Marshall, D.P., Popova, E., Treguier, A., Zanna, L., Yool, A., Yu, Y., Beadling, R., Bell, M.J., Kuhlbrodt, T., Arsouze, T., Bellucci, A., Castruccio, F., Gan, B., Putrasahan, D., Roberts, C.D., Van Roekel, L., Zhang, Q., 2020. Resolving and parameterising the ocean mesoscale in earth system models. Curr. Clim. Change Rep. 6, 137-152.

Hoegh-Guldberg, O.R., Cai, R., Poloczanska, E.S., Brewer, P.G., Sundby, S., Hilmi, K., Fabry, V.J., Jung, S., 2014. The ocean. In: Barros, V.R., Field, C.B., Dokken, D.J., Mastrandrea, M.D., Mach, K.J., Bilir, T.E., Chatterjee, M., Ebi, K.L., Estrada, Y.O., Genova, R.C., Girma, B., Kissel, A.N., MacCracken, S., Mastrandrea, P.R., White, L.
(Eds.), Climate Change 2014: Impacts, Adaptation, and Vulnerability. Part B: Regional Aspects. Contribution of Working Group II to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, pp. 1655-1731.

Holmes, R.M., Groeskamp, S., Stewart, K., McDougall, T.J., 2022. Sensitivity of a coarse-resolution global ocean model to spatially variable neutral diffusion. *J. Adv. Model. Earth Syst.* 14, e2021MS002914.

Jackson, L., Hughes, C.W., Williams, R.G., 2006. Topographic control of basin and channel flows: The role of the bottom pressure torques and friction. *J. Phys. Oceanogr.* 36, 1786-1805.

Jansen, M.F., Adcroft, A.J., Hallberg, R., Held, I.M., 2015a. Parameterization of eddy fluxes based on a mesoscale energy budget. *Ocean Model.* 92, 28-41.

Jansen, M.F., Adcroft, A., Khani, S., Kong, H., 2019. Toward an energetically consistent, resolution aware parameterization of ocean mesoscale eddies. *J. Adv. Model. Earth. Syst.* 1, 1-17.

Jansen, M.F., Held, I.M., Adcroft, A.J., Hallberg, R., 2015b. Energy budget-based backscatter in an eddy permitting primitive equation model. *Ocean Model.* 92, 15-26.

Jones, C.S., Abernathey, R.P., 2019. Isopycnal mixing controls deep ocean ventilation. *Geophys. Res. Lett.* 46, 13144-13151.

Juricke, S., Danilov, S., Koldunov, N., Oliver, M., Sein, D.V., Sidorenko, D., Wang, Q., 2020. A kinematic kinetic energy backscatter parametrization: From implementation to global ocean simulations. *J. Adv. Model. Earth Syst.* 12 (12), e2020MS002175.

Klocker, A., Marshall, D.P., 2014. Advection of baroclinic eddies by depth mean flow. *Geophys. Res. Lett.* 41, L60001.

LaCascé, J.H., Brink, K.H., 2000. Geostrophic turbulence over a slope. *J. Phys. Oceanogr.* 30 (6), 1305-1324.

Lee, M.-M., Coward, A.C., Nurser, A.G., 2002. Spurious diapycnal mixing of deep waters in an eddy-permitting global ocean model. *J. Phys. Oceanogr.* 32, 1522-1535.

Lévy, M., Iovino, D., Resplandy, L., Klein, P., Madec, G., Tréguier, A., Masson, S., Takahashi, K., 2012. Large-scale impacts of submesoscale dynamics on phytoplankton: Local and remote effects. *Ocean Modell.* 43-44, 77-93.

Lévy, M., Klein, P., Tréguier, A., Iovino, D., Madec, G., Masson, S., Takahashi, K., 2010. Modifications of gyre circulation by sub-mesoscale physics. *Ocean Model.* 34, 1-15.

Lindgren, F., Rue, H., Lindström, J., 2018. An explicit link between Gaussian fields and Gaussian Markov random fields: The stochastic partial differential equation approach. *J. R. Stat. Soc. Series B Stat. Methodol* 73, 423-498.

Madec, G., 2008. NEMO ocean engine. In: Note du Pôle de modélisation. (no. 27), Institut Pierre-Simon Laplace (IPSL).

Mak, J., Avdis, A., David, T.W., Lee, H.S., Na, Y., Yan, F.E., 2022a. On constraining the mesoscale eddy energy dissipation time-scale. *J. Adv. Model. Earth Syst.* 14, e2022MS003223.

Mak, J., Maddison, J.R., Marshall, D.P., Munday, D.R., 2018. Implementation of a geometrically informed and energetically constrained mesoscale eddy parameterization in an ocean circulation model. *J. Phys. Oceanogr.* 48, 2363-2382.

Mak, J., Maddison, J.R., Marshall, D.P., Ruan, X., Wang, Y., 2023. Scale-awareness in an eddy energy constrained mesoscale eddy parameterization. *J. Adv. Model. Earth. Syst.* 15, e2023MS003886.

Mak, J., Marshall, D.P., Madec, G., Maddison, J.R., 2022b. Acute sensitivity of global ocean circulation and heat content to eddy energy dissipation time-scale. *Geophys. Res. Lett.* 49, e2021GL097259.

Marshall, D.P., Maddison, J.R., Berloff, P.S., 2012. A framework for parameterizing eddy potential vorticity fluxes. *J. Phys. Oceanogr.* 42, 539-557.

Megann, A., 2018. Estimating the numerical diapycnal mixing in an eddy-permitting ocean model. *Ocean. Modell.* 121, 19-33.

Ni, Q., Zhai, X., Wang, G., Hughes, C.W., 2020a. Widespread mesoscale dipoles in the global ocean. *J. Geophys. Res. Oceans.* 125, e2020JC016479.

Ni, Q., Zhai, X., Wang, G., Marshall, D.P., 2020b. Random movement of mesoscale eddies in the global ocean. *J. Phys. Oceanogr.* 50, 2341-2357.

Nummelin, A., Isachsen, P.E., 2024. Parameterizing mesoscale eddy buoyancy transport over sloping topography. *ESS Open Archive*.

Perezhogin, P., Zanna, L., Fernandez-Granda, C., 2023. Generative data-driven approaches for stochastic subgrid parameterizations in an idealized ocean model. *J. Adv. Model. Earth Syst.* 15, e2023MS003681.

Pradal, M., Gnanadesikan, A., 2014. How does the Redi parameter for mesoscale mixing impact global climate in an Earth system model? *J. Adv. Model. Earth Syst.* 6, 586-601.

Redi, M.H., 1982. Oceanic isopycnal mixing by coordinate rotation. *J. Phys. Oceanogr.* 12, 1154-1158.

Roberts, M.J., Jackson, L.C., Roberts, C.D., Meccia, V., Docquier, D., Koenig, T., Ortega, P., Moreno-Chamarro, E., Bellucci, A., Coward, A.C., Drijfhout, S., Exarchou, E., Gutjahr, O., Hewitt, H., Iovino, D., Lohmann, K., Putrasahan, D., Schiemann, R., Seddon, J., Terray, L., Xu, X., Zhang, Q., Chang, P., Yeager, S.G., Castruccio, F.S., Zhang, S., Wu, L., 2020. Sensitivity of the Atlantic Meridional Overturning Circulation to model resolution in CMIP6 HighResMIP simulations and implications for future changes. *J. Adv. Model. Earth Syst.* 12, e2019MS002014.

Ruan, X., Couespel, D., Lévy, M., Mak, J., Wang, Y., 2023. Assessing the impact of ocean mesoscale eddy parameterizations on physical and biogeochemical responses: Diffusive closures in non-eddyding ocean models. *Ocean Model.* 183, 102204.

Séférían, R., Nabat, P., Michou, M., Saint-Martin, D., Voldoire, A., Colin, J., Decharme, B., Delire, C., Berthet, S., Chevallier, M., Sénési, S., Franchisteguy, L., Vial, J., Mallet, M., Joetzjer, E., Geoffroy, O., Guérémy, J., Moine, M., Msadek, R., Ribes, A., Rocher, M., Roehrig, R., Salas-y-Mélia, D., Sanchez, E., Terray, L., Valcke, S., Waldman, R., Aumont, O., Bopp, L., Deshayes, J., Éthé, C., Madec, G., 2019. Evaluation of CNRM Earth System Model, CNRM-ESM2-1: Role of Earth system processes in present-day and future climate. *J. Adv. Model. Earth. Syst* 11 (12), 4182-4227.

Smith, R.D., Gent, P.R., 2004. Anisotropic Gent-McWilliams parameterization for ocean models. *J. Phys. Oceanogr.* 34, 2541-2564.

Stewart, A.L., Thompson, A.F., 2013. Connecting Antarctic cross-slope exchange with Southern Ocean overturning. *J. Phys. Oceanogr.* 43 (7), 1453-1471.

Swearer, S.E., Treml, E.A., Shima, J.S., 2019. A Review of Biophysical Models of Marine Larval Dispersal. *CRC Press*.

Tokarska, K.B., Stolpe, M.B., Sippel, S., Fischer, E.M., Smith, C.J., Lehner, F., Knutti, R., 2020. Past warming trend constrains future warming in CMIP6 models. *Sci. Adv.* 6, eaaz9549.

Torres, R., Waldman, R., Mak, J., Séférían, R., 2023. Global estimation of the eddy kinetic energy dissipation from a diagnostic energy balance. *Geophys. Res. Lett.* 50, 2023GL104688.

Wang, Y., Stewart, A.L., 2018. Eddy dynamics over continental slopes under retrograde winds: Insights from a model inter-comparison. *Ocean Model.* 121, 1-18.

Waterman, S., Hoskins, B.J., 2013. Eddy shape, orientation, propagation, and mean flow feedback in Western boundary current jets. *J. Phys. Oceanogr.* 43, 1666-1690.

Waterman, S., Jayne, S.R., 2012. Eddy-driven recirculations from a localized transient forcing. *J. Phys. Oceanogr.* 42, 430-447.

Wei, H., Wang, Y., Mak, J., 2024. Parameterizing eddy buoyancy fluxes across prograde shelf/slope fronts using a slope-aware GEOMETRIC closure. *J. Phys. Oceanogr.* 54, 359-377.

Wei, H., Wang, Y., Stewart, A.L., Mak, J., 2022. Scalings for eddy buoyancy fluxes across prograde shelf/slope fronts. *J. Adv. Model. Earth. Syst.* 14, e2022MS003229.

Whittle, P., 1963. Stochastic processes in several dimensions. *Bull. Int. Stat. Inst* 40, 974-994.

Wilder, T., Munday, X.Z.D.R., Joshi, M., 2023. Constraining an eddy energy dissipation rate due to relative wind stress for use in energy budget-based eddy parameterisations. *Ocean Sci.* 19, 1669-1686.

Yankovsky, E., Bachman, S.D., Smith, K.S., Zanna, L., 2023. Vertical structure and energetic constraints for a backscatter parameterization of ocean mesoscale eddies. *ESS Open Archive*.