tmm4py: Global Ocean Biogeochemical Modeling in Python With the Transport Matrix Method

Abstract

Marine biogeochemical models are important tools in the quest to understand the cycling of chemical and biological tracers such as nutrients, carbon and oxygen, as well as key components of the Earth System Models used to project climate change. Historically, given the need for speed, global scale modeling has been performed in compiled languages like Fortran. However, as high level scripting languages such as Python and Julia gain popularity, the need for models and tools accessible from them has become imperative. This paper introduces tmm4py, a Python interface to a redesigned version of the Transport Matrix Method (TMM) software, a computationally efficient numerical scheme for “offline” simulation of marine geochemical and biogeochemical tracers. The TMM provides a convenient framework for developing and testing new biogeochemical parameterizations, as well as running existing complex models driven by circulations derived from state‐of‐the‐art physical models. tmm4py exposes all of the TMM library's functionality in Python, including transparent parallelization, allowing users to not only interactively use models written in compiled languages, but also develop complex models in pure Python with performance similar to compiled code. tmm4py enables users to exploit the large Python‐based scientific software ecosystem, including libraries for machine learning and deploying models on Graphics Processing Units. The various features of tmm4py are described and illustrated through practical examples, including a full‐fledged biogeochemical model written entirely in Python.

Introduction

Coupled physical-biogeochemical models of the ocean are an important part of the arsenal researchers use to understand the marine carbon cycle and its interaction with the climate system. Since they were first built many decades ago, these models have traditionally been implemented in compiled languages such as Fortran, as speed and the ability to exploit parallel architectures has always been critical to their practical use in a global setting. While the development of efficient “offline” approaches to biogeochemical modeling (e.g., DeVries & Primeau, 2011; Khatiwala, 2007; Khatiwala et al., 2005; Primeau, 2005) have simplified the task of writing and running global models, the language of choice has remained Fortran as both the language and compiler technology have evolved to keep up with the demanding requirements of increasing complexity and resolution. This is even as interpreted or scripting languages such as Python and Julia have come to the fore, completely displacing in many applications and fields compiled languages. This is somewhat ironic as scripting languages are only practically useful because the underlying libraries are implemented in C and Fortran (several decades ago in critical areas such as linear algebra). Nonetheless, the convenience and allure of high level scripting languages is undeniable and will only increase in the future.

This paper describes tmm4py, a new software aimed at bringing the convenience of Python to global ocean geochemical and biogeochemical modeling. tmm4py is based on the “Transport Matrix Method” (TMM), a type of “offline” ocean model in which passive tracers are advected and diffused using a precomputed circulation field from a general circulation model (GCM) (DeVries & Primeau, 2011; Khatiwala, 2007; Khatiwala et al., 2005; Primeau, 2005). The TMM represents tracer advection-diffusion as a sequence of sparse matrix-vector products, which can be performed much more quickly than running the full ocean GCM. The TMM and its implementation in C (Khatiwala, 2007, 2018) was originally envisioned as a way to “democratize” modeling by providing a convenient framework and software tools to rapidly develop and test new parameterizations or efficiently run existing models driven by circulation fields extracted from state-of-the-art GCMs. It has been used in numerous studies, ranging from constraining mean ocean temperature during the Last Glacial Maximum (Seltzer et al., 2024) to future changes in the biological carbon pump under global warming (Wilson et al., 2022), and from evaluating gas transfer parameterizations (Liang et al., 2013) to the first process model of neodymium isotopes (Siddall et al., 2008).

tmm4py provides a wrapper around the TMM software, redesigned along object oriented principles as a library, making it possible to write complex biogeochemical models in any language. By exposing all of the TMM's functionality, tmm4py combines the familiarity and interactivity of Python, while retaining the computational efficiency of the TMM. With tmm4py, models can be written entirely in Python, yet run at the speed of compiled code. In the following, the TMM and its redesigned implementation are briefly summarized, followed by a description of tmm4py and practical examples illustrating its features. While the focus of this article is on tmm4py, readers without an interest in writing Python code–existing TMM users in particular–may also find the description of the TMM library and the new capabilities enabled by tmm4py useful.

The Transport Matrix Method

In this section, the TMM is briefly described (see Khatiwala et al. (2005) and Khatiwala (2007) for details). The basic idea behind the TMM is that, when spatially discretized, the advective-diffusive transport of a tracer can be written as the product of a sparse matrix, that is, a matrix most of whose elements are zero, with a vector representing the tracer concentrations at the grid points of the underlying ocean model. The advantage of this formulation is that sparse matrix-vector products can be performed very efficiently on distributed memory computers and, increasingly, on Graphics Processing Units (GPUs; Bell and Garland (2008)). To accommodate the diversity of temporal discretization methods used in ocean GCMs, the TMM time-steps the following equation: 1 $${\mathbf{c}}^{n+1}={\mathbf{A}}_{i}^{n}\left({\mathbf{A}}_{e}^{n}{\mathbf{c}}^{n}+{\mathbf{q}}^{n}\right).$$ Here, $$\mathbf{c}$$ is the vector of tracer concentrations at the grid points of the GCM, that is, the 3-dimensional tracer field mapped to a vector; $$\mathbf{q}$$ represents sources and sinks, including surface and bottom (sediment) fluxes; $${\mathbf{A}}_{e}$$ the TM for the explicit-in-time component of advection-diffusion; $${\mathbf{A}}_{i}$$ the matrix for implicit transport (e.g., vertical mixing); and $$n$$ the time step index. For problems with multiple interacting tracers, one such equation is solved per tracer, with $$\mathbf{q}$$ in general a function of all tracers.

For problems in which tracers are subject to concentration boundary conditions (BC), the TMM can be formulated as: 2 $${\mathbf{c}}_{\mathrm{I}}^{n+1}={\mathbf{A}}_{i,\mathrm{I}}^{n}\left({\mathbf{A}}_{e,\mathrm{I}}^{n}{\mathbf{c}}_{\mathrm{I}}^{n}+{\mathbf{B}}_{e}^{n}{\mathbf{c}}_{\mathrm{B}}^{n}+{\mathbf{q}}_{\mathrm{I}}^{n}\right)+{\mathbf{B}}_{i}^{n}{\mathbf{c}}_{\mathrm{B}}^{n+1}$$ where, the subscripts “I” and “B” indicate interior and boundary values, respectively. The TMs that appear in this equation can be trivially constructed from the full TMs in Equation 1 by extracting appropriate rows and columns from the latter (Khatiwala, 2007).

It should be noted that there are other transport matrix-based approaches to solving the advection-diffusion equation (e.g., DeVries & Primeau, 2011; Primeau, 2005). These have generally been based on fixed or cyclostationary (repeating seasonal cycle) circulations. In contrast, the TMM allows for arbitrary time dependence, that is, the circulation can be fixed, cyclostationary or evolve over time as, for example, from a climate change simulation. Moreover, Equation 1 is applicable to and consistent with the wide majority of models and their numerics, as the TMs are extracted by “probing” the model with an automated, empirical procedure (Khatiwala et al., 2005), whereas alternative approaches are often hard-wired and tied to a specific ocean model.

The TMM Software

While the above equations can be time-stepped using any software, for efficiency the TMM has historically been implemented in C (Khatiwala, 2007, 2018) using PETSc (“Portable, Extensible Toolkit for Scientific Computation”; Balay et al., 1997, 2024), an open source suite of scalable solvers for linear and nonlinear equations. In addition to speed and a deep stack of software tools for managing matrices, vectors and other data structures in parallel, portability, long-term support (PETSc is now nearly 30 years old) and a forward-looking roadmap (GPU support has been steadily added) are some of the reasons for the decision to use PETSc to implement the TMM.

In the TMM, the forcing term $$\mathbf{q}$$ in Equations 1 and 2 can be either read from file ($${\mathbf{q}}_{\mathrm{f}}$$), imposed as a restoring term ($${\mathbf{q}}_{\mathrm{r}\mathrm{e}\mathrm{l}}=-\lambda \left(\mathbf{c}-{\mathbf{c}}_{\mathrm{r}\mathrm{e}\mathrm{l}}\right)$$), or computed ($${\mathbf{q}}_{\mathrm{e}\mathrm{f}}$$). Similarly, boundary conditions (Equation 2) can either be read from file or calculated. $${\mathbf{q}}_{\mathrm{f}}$$, $${\mathbf{q}}_{\mathrm{r}\mathrm{e}\mathrm{l}}$$ and BCs read from file are handled internally by the TMM, while $${\mathbf{q}}_{\mathrm{e}\mathrm{f}}$$ and calculated BCs require user-supplied “external forcing” functions. (This naming convention is inspired by MITgcm.) These C functions are in many cases simply wrappers around existing biogeochemical models written in Fortran. Amongst other models, wrappers exist for MITgcm (Dutkiewicz et al., 2005), BLING (Galbraith et al., 2010; Verdy & Mazloff, 2017), MOBI (Khatiwala et al., 2019; Schmittner & Somes, 2016)), MOPS (Kriest & Oschlies, 2015), MEDUSA (Yool et al., 2013) and PISCES (Aumont et al., 2015; Séférian et al., 2020), as well as for models simulating CFCs (Orr, Dutay, et al., 1999), noble gases (Nicholson et al., 2011, 2016; Seltzer et al., 2023), radiocarbon (Khatiwala et al., 2018; Orr, Najjar, et al., 1999), and various trace elements and their isotopes (e.g., de Souza et al., 2018; Siddall et al., 2008; Vance et al., 2017). These wrappers manage forcing data and map tracer and other fields between the TMM and model. Custom data structures and high level functions in the TMM software give users convenient ways to efficiently read and write in parallel 1-, 2-, and 3-dimensional fields common in biogeochemical modeling, for example, to read and interpolate in time seasonally-repeating surface wind speed. Other utility functions simplify the task of modeling biogeochemical processes that act in the vertical. All of this functionality can be used in serial or parallel without any knowledge of MPI programming.

The TMM software has a relatively “flat” structure composed of a main driver routine and a set of functions which time-step the above equations. Problem specific external forcing functions supplied by the user are called at the appropriate time during the time-stepping sequence. However, this design makes it difficult to couple to a language like Python. One possibility would have been to reimplement the TMM in Python using petsc4py (Dalcin et al., 2011), a module that provides deep integration with PETSc. However, this would have required rewriting a large amount of code developed and tested over many years with likely loss in efficiency, not to mention breaking compatibility with the many biogeochemical models already interfaced to the TMM. Instead, it was decided to redesign the existing software following PETSc's underlying object oriented philosophy (), while reusing much of the existing code. Briefly, PETSc solvers for different types of problems are encapsulated as classes that contain both data (the complete internal state of the algorithm) as well as pointers to functions that can operate on the data. The solvers can be applied to a specific problem through callback functions. For example, the SNES class solves nonlinear systems of equations $$\mathbf{F}(\mathbf{x})=\mathbf{0}$$. Users provide their own functions to compute $$\mathbf{F}(\mathbf{x})$$ and (optionally) its Jacobian matrix (or a function that computes the action of the Jacobian on a vector), and register them with a specific instance of the SNES object. PETSc's containerization technology makes it possible to attach to the object any type of problem-specific data required to evaluate the functions. An application code can contain multiple instances of the SNES class operating on different problems, each with their own associated data, with particular algorithms and algorithm parameters for each instance selectable at runtime. While conceptually simple, a sophisticated code infrastructure–arguably PETSc's “secret sauce”–is needed to enable this high level of encapsulation and abstraction.

Following these principles, a new solver class TMMState was created to encapsulate the full state of a biogeochemical model. (Like other PETSc classes, TMMState is a pointer to a C struct.) Users register their own functions to compute source terms and boundary conditions. Multiple instances of TMMState can exist within the same application code, each one corresponding to a difference instance of a biogeochemical model or an entirely different model altogether.

In the new TMM version, all relevant functions now take a TMMState object as an argument so that data encapsulation is maintained throughout, although user functions can choose to share data (e.g., forcing fields) between instances. Other classes were created to replace existing data structures that handle forcing data, timing and I/O. These include timer (StepTimer, PeriodicTimer, TimeDependentTimer), vector (PeriodicVec, TimeDependentVec) and array (PeriodicArray, TimeDependentArray) classes that provide high level interfaces for reading, writing and interpolating various types of data. For example, TimeDependentArray can be used to read 2-d forcing fields with arbitrary time dependence, while StepTimer deals with data discretely spaced (evenly or unevenly) in time (a typical use might be to accumulate and write out calendar month averages of a diagnostic).

With the new design the core TMM functionality can be separated from the driver. In particular, the former can be compiled into a static or dynamic library callable from programs written in C, Python and other languages. At the same time, despite the extensive refactoring, in all practical aspects the new code is almost completely backwards compatible in that existing users will encounter very little difference.

tmm4py: A Python Wrapper for the TMM

To create a bridge between the TMM and Python, a C extension module was built using the Cython compiler (Behnel et al., 2010, ). Cython allows developers to mix Python and C in the same source file (with extension .pyx), with any Python code being translated to optimized C code. The resulting .c file produced by Cython, which contains the code of a Python extension module, is then compiled by a C compiler into a dynamic library that can be imported directly into Python. Cython is used extensively in many Python projects (e.g., SciPy), and the petsc4py wrapper around PETSc is written in Cython.

Here, the main wrapper file is called tmm4py.pyx. It defines a Python class TMMState corresponding to the C TMMState class in the TMM library. To distinguish the two, the latter is hereafter called PetscTMMState. A schematic of the relationship between PETSc, TMM and tmm4py is shown in Figure 1. Only an essential subset of the data encapsulated in PetscTMMState is exposed to TMMState and whence in Python. (The variables exposed are defined in the corresponding tmm4py.pxd header file.) In the PetscTMMState C struct, fields corresponding to the tracer, source and boundary condition vectors are stored as PETSc Vec objects. These internally store the data as C arrays distributed across multiple processes, and also include (non-public) information about parallel layout, whether the data exist on CPU or GPU, and so forth. When a TMMState object is instantiated, tmm4py exposes those fields in Python as NumPy arrays whose memory is mapped to the corresponding underlying C arrays. If the code is run in parallel, only the “piece” of the array belonging to the current MPI process (rank) is made available to that process. In particular, if state is an instance of TMMState, then state.c is a list (of length number of tracers being simulated) of NumPy arrays containing the corresponding tracer values. Similarly, state.qef is a list of arrays for the user-computed source term. Table 1 lists the main variables in TMMState exposed by tmm4py.

[IMAGE OMITTED. SEE PDF]

Table 1 List of Variables Exposed by tmm4py in the TMMState Object

Variable name	Python data type	Description
c	list of NumPy arrays	Tracer values ($$\mathbf{c}$$ in Equation 1)
q	”	Source term ($$\mathbf{q}$$ in Equation 1 or $${\mathbf{q}}_{\mathrm{I}}$$ in Equation 2)
qef	”	User-computed “external forcing” source term ($${\mathbf{q}}_{\mathrm{e}\mathrm{f}}$$)
qf	”	Source term read from file ($${\mathbf{q}}_{\mathrm{f}}$$)
qrel	”	Restoring source term ($${\mathbf{q}}_{\mathrm{r}\mathrm{e}\mathrm{l}}$$)
cbc	”	Boundary values at the current time step ($${\mathbf{c}}_{\mathrm{B}}^{n}$$ in Equation 2)
cbf	”	Boundary values at the next time step ($${\mathbf{c}}_{\mathrm{B}}^{n+1}$$ in Equation 2)
config	dict	Dictionary of runtime options

The class also contains methods that can operate on TMMState objects. This includes functions for registering and evaluating user-specified functions to calculate external forcings and boundary conditions. These functions can be written in Python, C or Fortran, or a combination thereof. They can even be written in languages that interoperate with Python, although depending on how data are exchanged between the languages there may be some overhead. For example, models written in Julia can be used with tmm4py via modules such as JuliaCall ().

Since neither tmm4py nor the TMM library know anything at compilation time about the user function they will receive, to enable this flexibility callbacks are implemented as a two step process illustrated in Figure 2. First, the callback functions actually registered with the underlying PetscTMMState object by the TMM library are hardwired functions in tmm4py. It is the latter functions that are called by the TMM library during the time-stepping sequence and which in turn call the appropriate user functions. Second, in order for the hardwired functions in tmm4py to know which user function to call and with what arguments, during function registration, the user function, its arguments and the TMMState object itself are packed in a Python object and a pointer to it stashed by the TMM library in a PETSc container within the underlying PetscTMMState object. When a callback is triggered, this pointer is passed to the hardwired function, which unpacks the information and calls the user function with the TMMState object and other required arguments. Whew! This procedure may seem convoluted, but it ensures that tmm4py can handle arbitrary user functions and arguments. For external forcing functions written in Python, all variables required to compute source terms and boundary conditions are made available as Python data types (built-in, NumPy arrays, user defined classes, etc.). No additional code or intervention on the part of the user is required to communicate with the TMM library, and users can leverage existing Python modules, for example, PyCO2SYS (Humphreys et al., 2022) to solve the carbonate chemistry equations or gasex-python (D. Nicholson; ) to calculate gas transfer.

[IMAGE OMITTED. SEE PDF]

tmm4py.pyx also defines other classes corresponding to the TMM timer (StepTimer, PeriodicTimer, TimeDependentTimer), vector (PeriodicVec, TimeDependentVec) and array (PeriodicArray, TimeDependentArray) classes (see above). The underlying data are stored by the TMM in C arrays or PETSc Vecs and, as described further below, exposed to Python as NumPy arrays.

Note that for NumPy arrays mapped to underlying C arrays (directly or within PETSc objects), the memory is “owned” by the latter and garbage collection is performed when the arrays/objects are destroyed within the tmm4py module or TMM library. There is no deallocator method set for such NumPy arrays and thus deleting them will have no effect on the data stored in the C array. However, for certain classes (see below) a user can choose to use a NumPy array created within their Python code to store the underlying data. The constructor method for these classes take a NumPy array as an optional argument and point the underlying C array to the NumPy array's memory buffer. In such cases, memory is managed by the NumPy array and its deallocator method (called automatically when the array goes out of scope or when explicitly deleted) is set to free the C array pointer.

Examples

In the following, the various capabilities of tmm4py are demonstrated through a series of examples in which the model is implemented in Python. Given the practical nature of this paper, the section is organized as a tutorial with illustrative snippets of code. The complete code, run scripts, and detailed instructions for installing tmm4py and running the examples are provided as a supplement (see Open Research). Where relevant, Fortran versions of the models can also be run either from tmm4py or via the “classic” C TMM driver. They produce, within numerical roundoff, identical results.

Example 1: Age Tracer

We begin with a simple “hello world” example, the ideal age tracer (England, 1995; Thiele & Sarmiento, 1990), a passive tracer with a zero concentration boundary condition at the surface and an interior source term of “1” per unit time representing aging of water as it is transported by ocean circulation. In terms of Equation 2 this is modeled by setting $${\mathbf{c}}_{\mathrm{B}}=0$$ and $$\mathbf{q}$$ = 1. While this particular example doesn't actually require writing extra code–the TMM already has built-in capabilities for handing such cases (here, by reading from file the vector $$\mathbf{q}$$)–here we implement it as in a conventional ocean model via a user-calculated external forcing term.

External forcing in the TMM is enabled with the -external_forcing runtime flag, which instructs the TMM to call the registered user callback functions during the time-stepping sequence. These can be standalone functions but it is convenient to encapsulate them in a Python class as shown in Listing 1 for the age model. The listing shows three functions, to initialize, calculate and finalize, respectively, the forcing. (Specifying these is mandatory, even if they don't perform any action.)

1

Listing

Implementation of a Python class for the ideal age tracer.

class Age:  # Constructor  def __init__(self,DeltaT=28800.0):    self.dt=DeltaT  # Initialize the model  def iniAgeTracer(self,tc,Iter,state,*args,**kwargs):    pass  # Main source computation function  def calcAgeTracer(self,tc,Iter,iLoop,state,*args,**kwargs):   # time−discrete age tracer source term ([second]*[year/second])   state.qef [0][:] = self.dt*1.0/(86400.*365.)  # Clean up as necessary  def finalizeAgeTracer(self,tc,Iter,state,*args,**kwargs):    pass

As mentioned above, tmm4py automatically makes available the underlying tracer and source vectors as NumPy arrays. In the line state.qef[0][:]=self.dt*1.0/(86400.*365.), qef is a list of length number of simulated tracers (1 here) of NumPy arrays corresponding to the “piece” of the external forcing vector $${\mathbf{q}}_{\mathrm{e}\mathrm{f}}$$ belonging to this MPI rank. state.qef[0][:]=… thus sets all elements of the first source array (index 0) to the forcing term.

The Age module can then be used in a driver script that time-steps the TMM tracer equations. A typical script would carry out the following steps:

• Import required modules

• Initialize the PETSc and TMM libraries

• Create a biogeochemical model instance

• Instantiate a TMMState object and set various runtime options corresponding to that state

• Register model functions with the state object

• Time-step the tracer equations

These steps are common to all biogeochemical models and a driver script only need be written once; we simply import the relevant model module and register its functions with a TMMState object.

A fully functional driver script is shown in Listing 2. It can be run without modification either in serial:

python driver.py

or in parallel, for example:

mpiexec −np 48 python driver.py

2

Listing

Driver script to simulate the Age model.

# Import required modules including the tmm4py wrapper import petsc4py import sys import tmm4py as TMM from Age import Age # Initialize the PETSc and TMM libraries petsc4py.init(sys.argv) Iter0,maxSteps,time0,deltaTClock=TMM.initialize() # Create an instance of the Age model model = Age() # Instantiate a TMMState object and set various runtime options state = TMM.TMMState() state.create() state.setFromOptions(doOutput=True) # Register Age model callback functions with the TMMState instance state.setIniExternalForcingFunction(model.iniAgeTracer) state.setCalcExternalForcingFunction(model.calcAgeTracer) state.setFinExternalForcingFunction(model.finalizeAgeTracer) # Time−step tracer equations for iLoop in range(1,maxSteps + 1):  tc = time0 + deltaTClock*(iLoop − 1)  tf = time0 + deltaTClock*iLoop  Iterc = Iter0 + iLoop − 1  TMM.updateTMs(tc) # update transport matrices  state.forcingUpdate(tc,Iterc,iLoop) # update forcing and BCs  state.timeStep(tc,Iterc,iLoop) # step forward tracers  state.timeStepPost(tf,Iterc,iLoop) # post time step functions  state.output(tf,Iterc,iLoop) # write output

While not especially useful in this example, it is also possible to instantiate and simulate multiple (independent) copies of the model, each corresponding to a different instance of TMMState. To facilitate this, we exploit PETSc's extensive set of functions for processing runtime options (exposed in Python by petsc4py). At runtime, PETSc generates a database of options passed via the command line, read from a file, specified in the PETSC_OPTIONS environment variable, or embedded within the code itself. This database can then be queried within the code to set parameter values. For instance, in the TMM library the line:

PetscOptionsGetStringArray(NULL,prefix,”−i”,state−>iniFile, &maxValsToRead,&flg)

queries the database for the runtime option “-i” and corresponding list of initial condition file names (specified as, e.g., -i file1,file2). To read different sets of initial condition files for different TMMState instances, it is not necessary to modify the code. Conveniently, PETSc gives users the ability to tag options with a “prefix” (passed as an argument to database query functions, as in the above line of code). Thus, initial conditions for a state with prefix “state1_” can be specified on the command line with, for example, -state1_i file1,file2.

Returning to the age example, the steps involved are:

• Access the PETSc options database

• Query the database for a list of prefixes that tag the different TMMState instances (specified at runtime with, for example, -prefixes age1_,age2_)

• Create the requested number of copies of Age

• Instantiate TMMState objects and register functions for each TMMState/Age pair

Listing 3 shows a modified version of the previous driver code implementing these steps.

3

Listing

Modified driver to run multiple instances of Age.

# Import required modules including the tmm4py wrapper import petsc4py import sys import tmm4py as TMM from Age import Age from petsc4py import PETSc # Initialize the PETSc and TMM libraries petsc4py.init(sys.argv) Iter0,maxSteps,time0,deltaTClock=TMM.initialize() OptDB=PETSc.Options() # Get options database as a dictionary if OptDB.hasName(”prefixes”):    prefixes=OptDB.getString(”−prefixes”).split(”,”) else:    prefixes=[””] # Create desired number of instances of Age model n=len(prefixes) models]= [Age() for _ in range(n) # Instantiate corresponding TMMState objects states=[None]*n for i in range(n):    states[i]=TMM.TMMState()    states[i].create()    states[i].setFromOptions(prefix=prefixes[i],doOutput=True)    states[i].setIniExternalForcingFunction(models[i].iniAgeTracer)    states[i].setCalcExternalForcingFunction(models[i].calcAgeTracer)    states[i].setFinExternalForcingFunction(models[i].finalizeAgeTracer) # Time−step tracer equations for iLoop in range(1,maxSteps + 1):    tc=time0 + deltaTClock * (iLoop − 1)    tf=time0 + deltaTClock * iLoop    Iterc=Iter0 + iLoop − 1   TMM.updateTMs(tc) # update transport matrices    for i in range(n):        states[i].forcingUpdate(tc,Iterc,iLoop) # update forcing and BCs        states[i].timeStep(tc,Iterc,iLoop) # step forward tracers        states[i].timeStepPost(tf,Iterc,iLoop) # post time step functions        states[i].output(tf,Iterc,iLoop) # write output

Incidentally, all the registration functions (setIniExternalForcingFunction, etc) take optional arguments that are passed to the corresponding callback function when it is invoked. Thus, if a function is registered with:

state.setIniExternalForcingFunction(model.iniAgeTracer,3,’a’,p1=17.0,p2=3.0),

as with any Python function, the first two optional arguments are packed into a tuple (3, ’a’) and the next two keyword arguments into a dictionary {’p1’: 17.0, ’p2’: 3.0} and passed to args and kwargs, respectively, in iniAgeTracer (see Listing 1). This can be useful, for example, for setting model parameter or running multiple model instances with different parameter values (the optional arguments are bound to the TMMState instance registering the function; see Figure 2).

Example 2: Protactinium and Thorium Isotopes

We next consider a slightly more complex model involving the radioactive tracers ²³¹Pa and ²³⁰Th. These tracers are produced by uranium decay and scavenged onto sinking biogenic and lithogenic particles (Henderson & Anderson, 2003). The ratio of ²³¹Pa to ²³⁰Th is commonly used as a paleoproxy for the strength of the Atlantic Meridional Overturning Circulation (Henderson & Anderson, 2003; McManus et al., 2004; Yu et al., 1996). The tracers are widely implemented in ocean and climate models used for paleoclimate studies (e.g., Gu & Liu, 2017; Missiaen et al., 2020; Rempfer et al., 2017; Sasaki et al., 2022; van Hulten et al., 2018). Pa/Th models share a number of features with more complex biogeochemical models, in particular the fact that both involve external forcing data and processes that act in the vertical, making this a good example to illustrate how tmm4py can be used to implement such models. In the following, only snippets of code relevant to demonstrating these capabilities are shown. Readers can refer to the provided code for the full listing.

For this example, the (Fortran) Pa/Th module in MOBI (Khatiwala et al., 2019; Schmittner & Somes, 2016) was translated into Python and implemented as a class. A full description of this MOBI module is currently being written for submission. Briefly, it is based on Siddall et al. (2005), and treats adsorption onto and desorption from 4 different particle types (particulate organic carbon, opal, calcium carbonate and dust) via the reversible scavenging parameterization of Bacon and Anderson (1982).

As with the previous age tracer example, the PaTh class implements methods for model initialization, calculation of the source term, and deleting objects at the end of the integration. As is typical of biogeochemical models, the initialization function reads grid and forcing data. Simple binary files can be read using built-in Python functions, for example:

# Thickness of grid boxes

with open(”drF.bin”,’rb’) as f:

nzmax=np.fromfile(f,dtype=np.dtype(’>i4’),count=1)[0]

drF=np.fromfile(f,dtype=np.dtype(’>f8’),count=p.nzmax)

Two- and three-dimensional fields can be read using functions such as tmm4py.loadVecIntoArray, which ensures that data are distributed with the appropriate parallel layout for that field and returns the “piece” belonging to the MPI rank as a NumPy array:

# Read 3−d particle concentration field

localPaTh_pom=tmm4py.loadVecIntoArray(”PaTh_pom.petsc”)

(For code readability, variables whose names start with “local” are the local piece of arrays distributed over multiple processes.)

Seasonally-varying forcing fields can be handled by a combination of a PeriodicTimer object and PeriodicArray (2-d) or PeriodicVec (3-d) objects, which read bracketing slices from file and (internally) store them in (non-public) C arrays or PETSc Vecs. Data with arbitrary time dependence can be handled with the corresponding TimeDependentTimer, TimeDependentArray and TimeDependentVec classes. For example, for a seasonally-varying surface field (dust flux here) we first create a PeriodicTimer object:

biogeochemTimer=tmm4py.PeriodicTimer()

biogeochemTimer.create(prefix=”periodic_biogeochem_”)

This is followed by creation of a PeriodicArray object, which in this example uses a user-supplied NumPy array as the memory buffer into which to place the data interpolated in time:

# Memory buffer of length=number of surface points on this MPI rank

localdustdep=np.zeros(lNumProfiles)

localdustdepp=tmm4py.PeriodicArray()

localdustdepp.create(lNumProfiles,buf=localdustdep)

Supplying a NumPy array is more for code readability than out of necessity as the interpolated data can also be directly accessed from the PeriodicArray.arr NumPy array. Note that timer objects only keep track of time; the same object can be shared between fields with the same temporal variation.

We then call the PeriodicArray.interp method to interpolate the field in time:

# Interpolate the dust flux to time t

localdustdepp.interp(t,biogeochemTimer,”dust_dep_”)

The interp function uses information in the timer object to read bracketing slices from file(s) as needed, before interpolating to the desired time. Here, the results are available in both localdustdepp.arr and the user-supplied localdustdep NumPy array (they share the same memory).

The computation of the source term $${\mathbf{q}}_{\mathrm{e}\mathrm{f}}$$ involves looping over vertical profiles (columns). When run with the -use_profiles option, the TMM library assumes that the tracer vectors are organized as a sequence of vertical profiles “stacked” one after the other. (It is up to the user to ensure that transport matrices and any forcing data are similarly arranged.) It then partitions the data across MPI processes ensuring that profiles aren't “broken” across processes. To access a vertical profile on a given MPI rank we need to know the local starting index of that profile and its length. Once these are known, we can then loop over each profile, “pulling out” as necessary the piece of the local tracer, source and forcing arrays corresponding to that profile and performing calculations with them. This pattern is common enough in biogeochemical models that to facilitate it, tmm4py exposes all the relevant parameters (computed within the TMM library) as a dictionary that can be accessed by calling tmm4py.getProfileConfig(). Listing 4 shows how this function is used in the PaTh module to perform a loop over columns.

4

Listing

Code showing how to access and use layout information when tracer and other vectors are arranged as a sequence of stacked vertical profiles (columns).

# Fetch dictionary of profile parameters profileConfig=tmm4py.getProfileConfig() # Profile parameters for this MPI rank # Total length of tracer vector lSize=profileConfig[’lSize’] # Number of profiles lNumProfiles]=profileConfig[’lNumProfiles’] # NumPy array of profile lengths lProfileLength=profileConfig[’lProfileLength’] # NumPy array of starting indices lStartIndices=profileConfig[’lStartIndices’] # Loop over profiles for ip in range(lNumProfiles):      nzloc=lProfileLength[ip] # length of profile      kl=lStartIndices[ip] # starting index      # NumPy arrays pointing to data for this vertical profile      # Tracers      Pa=state.c[0][kl:kl+nzloc]      Th=state.c[1][kl:kl+nzloc]      # Source terms      srcPa=state.qef[0][kl:kl+nzloc]      srcTh=state.qef[1][kl:kl+nzloc]      # Forcing      POM=localPaTh_pom[kl:kl+nzloc]      #etc      # Pass arrays for this profile to the main computational function      FPa_bot,FTh_bot=protac_thor_driver(nzloc,…,Pa,Th,POM,…)

Lastly, after time-stepping is complete, in the finalize function we delete any objects created by the model. Here, we call the destroy method on the PeriodicArray localdustp which releases memory (within the TMM library) for internal arrays used to store bracketing slices, temporary work space, and so on:

localdustp.destroy()

It is not necessary to call the finalize function explicitly; it is automatically invoked by the TMM library at the end of the program, when it also frees up memory internally allocated by PETSc.

Example 3: MOPS Biogeochemical Model

As a final example to show how tmm4py can be used to implement more complex biogeochemical processes, the MOPS model (Kriest & Oschlies, 2015) was translated from Fortran to Python. (The conversion, which took several therapeutic hours, was performed “by hand” without the aid of AI tools such as Copilot.) MOPS simulates the cycling of 9 biogeochemical tracers within the water column, including phytoplankton, zooplankton, and detritus; dissolved inorganic phosphate, nitrate and carbon; dissolved oxygen; alkalinity; and dissolved organic phosphate and nitrate. The version used here–“IronMOPS”–additionally simulates dissolved and particulate iron (I. Kriest, pers. comm.). It should be emphasized that the point here and in the previous example is not that existing models can or should be translated, but to demonstrate that such models can be easily implemented in a language like Python with the tools provided by tmm4py.

As is common in such models, MOPS has multiple subroutines spread over several source files that perform different tasks, such as solving the carbonate chemistry equations, calculating air-sea fluxes of CO₂, and so forth. In the Python version, this modular structure was largely maintained, based around a MOPS class similar to the PaTh class above. MOPS is obviously more complex, requiring many more forcing fields of different types (e.g., wind speed, iron flux from dust and hydrothermal vents, 3-d temperature and salinity, etc., all varying seasonally). But it is no more difficult in practice to implement.

MOPS also requires additional grid information such as the surface area of grid boxes for computing air-sea CO₂ fluxes and river runoff:

# Read 2−d surface data and return piece corresponding

# to this MPI rank as a NumPy array

localdA=tmm4py.readProfileScalarData(’dA.bin’)

Many models (MOPS included) additionally require the total surface area of the ocean. Since localdA above only contains the grid box areas for the current MPI process, a simple sum over the elements of localdA is not sufficient; we must also sum over all MPI processes. Users have the option of implementing such operations with the mpi4py Python module, for example:

# Load mpi4py module

from mpi4py import MPI

comm=MPI.COMM_WORLD

# First get total area for this MPI rank …

localA=np.sum(localdA)

# … and then sum over all MPI ranks to get the global area

totalA=comm.allreduce(localA,op=MPI.SUM)

Alternatively, they can use one of several convenience functions tmm4py provides for such common situations:

totalA=tmm4py.globalSum(localdA)

As another example, consider a biogeochemical model coupled to a well-mixed atmospheric box driven with CO₂ emissions. To prognostically time-step the atmosphere the globally-integrated air-sea CO₂ flux is needed. This can be calculated in two ways. Suppose localco2airseaflux and localdA are NumPy arrays of the air-sea CO₂ flux and surface area, respectively, on a MPI process. We can either use MPI:

# First sum the flux on this MPI rank …

localFocean=np.dot(localco2airseaflux,localdA)

# … and then sum over MPI ranks to get global air−sea CO2 flux

Focean=comm.allreduce(localFocean,op=MPI.SUM)

Or, we let tmm4py do it:

Focean=tmm4py.dotProd(localco2airseaflux,localdA)

The function tmm4py.dotProd is useful in other contexts too, for instance to calculate the global inventory of tracers in the model:

# Load array of grid box volumes

vol=tmm4py.loadVecIntoArray(”volumes.petsc”)

# Number of tracers for this TMMState instance

numTracers=state.config[’numTracers’]

# Compute global inventory of each tracer and return as a list

I=[tmm4py.dotProd(state.c[itr],vol) for itr in range(numTracers)]

In general the tmm4py functions will be more efficient as they use compiled code and optimized low level libraries.

MOPS's main computational routine BGC_MODEL is similar to the one for Pa/Th and other biogeochemical models in that it operates on a single vertical profile. Also, as is typical, this routine internally takes several small time steps to apply the biogeochemical source terms it computes to evolve the tracer field over one tracer advection-diffusion time step. That is, BGC_MODEL does not set the source arrays in state.qef, but directly modifies the tracer arrays in state.c. The loop over profiles in MOPS's external forcing function is more or less identical to that for Pa/Th (Listing 4) and not shown here for brevity.

Like any biogeochemical model, MOPS also computes diagnostics such as net primary production and respiration. Since the calculations are local to each MPI process, they can be easily handled by NumPy arrays and standard Python constructs. However since the arrays represent global distributed fields, to save them to file we call one of several functions in the TMM library exposed by tmm4py to perform parallel I/O. In particular, a NumPy array representing a distributed field (2-d or 3-d) can be written with a single function call:

tmm4py.writeArrayToVec(array,filename,…).,

As discussed above in the age tracer example (see Listing 3), the TMM can handle multiple instance of TMMState, each of them associated with a different model. For instance, to simultaneously simulate MOPS and Pa/Th we can do:

from MOPS import MOPS

from PaTh import PaTh

prefixes=[”mops_”,”path_”]

numStates=len(prefixes)

models=[MOPS(),PaTh()]

The prefixes can then be used to specify runtime options for each model by prefacing TMM options with the corresponding prefix. For example, to specify output file names (“-o”):

−mops_o po4.petsc,dop.petsc,oxy.petsc, …

−path_o Pa.petsc,Th.petsc

Interfacing tmm4py to Models Written in Fortran

The examples above show how models can be written purely in Python. Many users though will want to use existing models written in Fortran. Such models communicate with the TMM via an intermediate C wrapper (see Figure 1). As mentioned above, wrappers exist for a large number of models and these can continue to be used as before, either with the “classic” TMM C driver or from Python with tmm4py. While an additional communication layer is required to enable the latter, it allows the use of other Python libraries and a degree of interactivity not possible with the C version. (The same runtime options can be used whether a model is run from C or Python since those are handled by the back-end TMM/PETSc libraries; thus users considering switching to tmm4py won't need to modify run scripts written for the C driver.)

The general steps for using a Fortran model from tmm4py are as follows:

• Compile the (Fortran) model and (C) wrapper into a dynamic library.

• Write Cython “gateway” functions accessible from Python to call the wrapper functions. A typical implementation of such an interface is shown below:

  from petsc4py.PETSc import Error

  from tmm4py cimport PetscTMMState,TMMState

  # Get the wrapper function signature

  cdef extern from “tmm_external_forcing.h”:

  PetscErrorCode calcExternalForcing(…)

  # Python gateway function to C wrapper function

  def calcExternalForcingFn(tc,Iter,iLoop,TMMState State):

  # Call the wrapper function

  ierr=calcExternalForcing(…);

• Compile the above interface with Cython and compile the resulting .c file together with the model dynamic library into a Python extension module. This module can be directly imported into Python and referenced in exactly the same way as the previous pure Python examples. Thus, if the module is imported as model, the gateway function shown above can be registered with:

  state.setCalcExternalForcingFunction(model.calcExternalForcingFn)

The compilation steps can be easily automated and since all wrapper functions have the same signatures regardless of the underlying model, the gateway interface only needs to be written once. model_tmm_interface.pyx supplied with tmm4py contains a complete set of gateway functions that will suffice for most users. The above steps are implemented in the provided examples and demonstrate the use of the Fortran versions of Pa/Th and MOPS from tmm4py.

Performance of Biogeochemical Models Written in Python

While the convenience of writing complex models in a familiar language such as Python is appealing, it naturally comes at the cost of performance. Code written in Fortran and other compiled languages can run significantly faster; much of the core scientific Python stack, from NumPy to SciPy, is written in C or depends on libraries like BLAS and LAPACK written in C or Fortran. Scripting languages are generally only competitive when code can be vectorized to exploit low level compiled computational kernels. Unfortunately, biogeochemical models routinely involve multiple nested loops over vertical profiles that access array elements and are thus difficult to vectorize. Indeed, a comparison of the Python and Fortran versions of MOPS reveals that the former is about 250 times slower, making it practically unusable. (While more seasoned Python programmers than this author will no doubt be able to eke out better performance, it is highly unlikely to come anywhere close to Fortran.)

Fortunately, software tools such as “just-in-time” (JIT) compilers can translate interpreted code at runtime to machine code that executes much faster. Numba (Lam et al., 2015, ), for example, is a JIT compiler for Python especially built for code involving NumPy arrays and loops. Numba provides a set of “decorators” that can be applied to user functions. Numba works best on functions with clear, relatively simple inputs and outputs and data types it understands. Functions that take complex or custom objects as arguments are likely to be rejected by the Numba compiler. Thus, it cannot be applied to the external forcing functions directly. Nor is it necessary to, since most code in those functions involve calls to compiled functions (e.g., to read or interpolate fields). The exception is calls to core biogeochemical routines, for example, BGC_MODEL in MOPS. As long as these functions have a clean interface primarily involving NumPy arrays and basic Python data types, they can benefit tremendously from JIT compilation.

Applying Numba is as simple as:

# Import the decorator

from numba import jit

# Decorate the function to accelerate

@jit(nopython=True)

def BGC_MODEL(…):

# rest of function …

Other functions that can be similarly decorated and accelerated include the carbonate chemistry solver. The Numba-decorated version of MOPS is $${\sim} $$100 times faster than the original version, making it only between 2 and 3 times slower than the Fortran version.

Arguably, as implemented, the code does not optimally exploit JIT compilation. The main computational routine BGC_MODEL operates on a single column, which is a direct translation of the Fortran version. Each time a JIT-compiled function is called with different arguments there is some overhead that can add up. A simple refactoring to move the loop over columns to within BGC_MODEL speeds-up the JIT-compiled version by a further factor of 3, essentially eliminating the gap between the Python and Fortran versions. This experience illustrates that to make optimal use of libraries like Numba some redesign of the algorithm may be necessary.

It should be noted that other Python modules such as JAX (Bradbury et al., 2018) also have similar capabilities. Numba was chosen to illustrate JIT compilation since it operates on NumPy arrays directly by accessing their memory buffer, making it especially easy to apply and work with many existing Python modules. On the other hand, libraries like JAX don't directly work on NumPy arrays but provide their own array object with a “NumPy-inspired” interface. Moreover, in JAX's case, these arrays are immutable, don't generally support operations like slicing and assignment, and must be manipulated via JAX's API. However, for applications requiring differentiable code (e.g., for machine learning) that run on GPUs, JAX may be especially effective. Ultimately, the choice of a particular library will depend on the context and users' preferences.

Summary and Future Directions

This paper describes tmm4py, a new software to enable global scale marine geochemical and biogeochemical modeling in Python. It is built on top of the TMM software, which has been redesigned as an object oriented library that can be used from other programs and languages. The tmm4py wrapper exposes the TMM's extensive functionality, including transparent parallelization, allowing users to write complex models entirely in Python while exploiting the deep pool of software developed for that language. Models written in Fortran can continue to be used as before through tmm4py but with added benefits over the “classic” C TMM version. One of these is the ability to run biogeochemical models interactively. Unlike its C version or conventional models, the Python driver script (e.g., Listing 2) can be executed a line at a time in an interactive Python session, that is, one can take a single time step, examine or plot the solution, take another step, and so on. Users can even directly manipulate model fields–whether written in Python or Fortran–or call model functions individually. These capabilities may be especially useful when debugging a model or as a pedagogical tool.

With tmm4py it is also possible to build more complex algorithms and workflows. For example, with a few dozen lines of code it can be used to perform sensitivity analysis of biogeochemical models written in either Python or Fortran by leveraging existing libraries such as JAX for automatic differentiation of Python functions. Another example is the efficient spin-up of seasonally-forced biogeochemical models using recently proposed algorithms such as Anderson Acceleration (Khatiwala, 2023, 2024). This method treats the model as a “black box”, typically communicating with it via restart files. However, because the scheme is implemented in Python and parallelized via MPI, it can share the same processor layout and directly interact with tmm4py without the overhead of reading/writing files. In fact, it is already planned to incorporate the functionality to spin-up models into tmm4py that can be enabled via a simple runtime switch.

tmm4py also makes it possible to run multiple models or instances of the same model simultaneously. As an example of where this might be useful, consider a “coupled” simulation in which particle sinking speeds and concentration fields computed by a biogeochemical model are passed to the Pa/Th module (instead of reading them from file). Another example is running multiple instances of the same model with different parameter values for calibration or sensitivity analysis. The advantage here is that the cost of reading/interpolating transport matrices and forcing fields is amortized over several model instances.

Global biogeochemical models have historically been implemented in Fortran to maximize performance. However, with just a couple of lines of extra code, just-in-time compilers can dramatically accelerate models written in Python, significantly reducing or even completely eliminating the speed advantage of models written in compiled languages. Such compilers can even target GPUs, which opens up the possibility of running models on such devices. Numba, for instance, directly compiles a subset of Python code into CUDA kernels that execute on NVIDIA devices (with an appropriate decorator a function can be JIT compiled for use on both CPU and GPU, thus allowing the same code to be reused). With compiled languages, the GPU software stack is largely skewed toward C and C++, making it much more challenging to port Fortran models to run on GPUs. (PETSc's GPU support is well developed (Mills et al., 2021), and the TMM library can already exploit it for aspects such as sparse matrix-vector multiplication.) JAX, CuPy (Okuta et al., 2017, ) and cuNumerica (Bauer & Garland, 2019, ) are other Python libraries geared toward making it easier to deploy code on GPUs.

Not that the task is completely trivial. Each library has a different API, which limits portability and risks lock-in. Languages like CUDA impose many constraints on what code running on the GPU can do. For example, it is not possible to create arrays within functions that run on devices; all arrays–even those used to hold data temporarily (as is typical in models)–must be created on the host and passed as arguments to the kernel function (which can also be challenging as CUDA limits the number of arguments that can be passed). Memory management thus requires careful consideration to minimize costly data transfer between the CPU and GPU, as does algorithm design. Many of the aforementioned libraries are more likely to be useful on problems with vectorized array operations or simple element wise operations on arrays for which users can write kernel functions for execution on GPUs. Computing complex biogeochemical interactions that involve, e.g., multiple array elements, may not easily fit into that paradigm. Nevertheless, these tools do lower the barrier to entry, and tmm4py provides a potential avenue for using them to perform biogeochemical modeling on this new architecture.

Lastly, while the focus of this paper has been on a Python interface, the design of the TMM library makes it suitable for building wrappers in other languages (e.g., Julia).

Acknowledgments

I am grateful to Jamie Carr for generously sharing his Python expertise, and to Iris Kriest for making available the source code of the latest version of her MOPS model. Thank you, too, to the two reviewers for their constructive comments that greatly improved the manuscript. I would like to acknowledge high-performance computing support from the Derecho system () provided by the NSF National Center for Atmospheric Research (NCAR), sponsored by the National Science Foundation.

Data Availability Statement

The TMM, tmm4py and example code, as well as all data needed to run the examples are archived in Khatiwala (2025) (). Also contained therein are detailed instructions for installing the software and required libraries, and for running the examples. The latest version of the code is available from .

References

Aumont, O., Ethé, C., Tagliabue, A., Bopp, L., & Gehlen, M. (2015). PISCES‐v2: An ocean biogeochemical model for carbon and ecosystem studies. Geoscientific Model Development, 8(8), 2465–2513. https://doi.org/10.5194/gmd‐8‐2465‐2015

Bacon, M. P., & Anderson, R. F. (1982). Distribution of thorium isotopes between dissolved and particulate forms in the deep sea. Journal of Geophysical Research, 87(2045), 216–237. https://doi.org/10.1029/JC087iC03p02045

Balay, S., Abhyankar, S., Adams, M. F., Benson, S., Brown, J., Brune, P., et al. (2024). PETSc/TAO users manual (Technical Report No. ANL‐21/39 ‐ Revision 3.22). Argonne National Laboratory. https://doi.org/10.2172/2205494

Balay, S., Gropp, W. D., McInnes, L. C., & Smith, B. F. (1997). Efficient management of parallelism in object oriented numerical software libraries. In E. Arge, A. M. Bruaset, & H. P. Langtangen (Eds.), Modern software tools in scientific computing (pp. 163–202). Birkhauser Press.

Bauer, M., & Garland, M. (2019). Legate NumPy: Accelerated and distributed array computing. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (Denver, Colorado) (SC ’19). Association for Computing Machinery. https://doi.org/10.1145/3295500.3356175

Behnel, S., Bradshaw, R., Citro, C., Dalcin, L., Seljebotn, D. S., & Smith, K. (2010). Cython: The best of both worlds. Computing in Science & Engineering, 13(2), 31–39. https://doi.org/10.1109/MCSE.2010.118

Bell, N., & Garland, M. (2008). Efficient sparse matrix‐vector multiplication on CUDA (Technical Report). NVIDIA Technical Report NVR‐2008‐004.

Bradbury, J., Frostig, R., Hawkins, P., Johnson, M. J., Leary, C., Maclaurin, D., et al. (2018). JAX: Composable transformations of Python+NumPy programs. Retrieved from http://github.com/jax‐ml/jax

Dalcin, L. D., Paz, R. R., Kler, P. A., & Cosimo, A. (2011). Parallel distributed computing using Python. Advances in Water Resources, 34(9), 1124–1139. https://doi.org/10.1016/j.advwatres.2011.04.013

de Souza, G. F., Khatiwala, S., Hain, M. P., Little, S. H., & Vance, D. (2018). On the origin of the marine zinc–silicon correlation. Earth and Planetary Science Letters, 492, 22–34. https://doi.org/10.1016/j.epsl.2018.03.050

DeVries, T., & Primeau, F. (2011). Dynamically and observationally constrained estimates of water‐mass distributions and ages in the global ocean. Journal of Physical Oceanography, 41(12), 2381–2401. https://doi.org/10.1175/JPO‐D‐10‐05011.1

Dutkiewicz, S., Follows, M., & Parekh, P. (2005). Interactions of the iron and phosphorus cycles: A three‐dimensional model study. Global Biogeochemical Cycles, 19(1), GB1021. https://doi.org/10.1029/2004GB002342

England, M. H. (1995). The age of water and ventilation timescales in a global ocean model. Journal of Physical Oceanography, 25(11), 2756–2777. https://doi.org/10.1175/1520‐0485(1995)025<2756:taowav>2.0.co;2

Galbraith, E. D., Gnanadesikan, A., Dunne, J. P., & Hiscock, M. R. (2010). Regional impacts of iron‐light colimitation in a global biogeochemical model. Biogeosciences, 7(3), 1043–1064. https://doi.org/10.5194/bg‐7‐1043‐2010

Gu, S., & Liu, Z. (2017). 231Pa and 230Th in the ocean model of the Community Earth System Model (CESM1.3). Geoscientific Model Development, 10(12), 4723–4742. https://doi.org/10.5194/gmd‐10‐4723‐2017

Henderson, G. M., & Anderson, R. F. (2003). The U‐series toolbox for paleoceanography. Reviews in Mineralogy and Geochemistry, 52, 493–531. https://doi.org/10.2113/0520493

Humphreys, M. P., Lewis, E. R., Sharp, J. D., & Pierrot, D. (2022). PyCO2SYS v1.8: Marine carbonate system calculations in Python. Geoscientific Model Development, 15(1), 15–43. https://doi.org/10.5194/gmd‐15‐15‐2022

Khatiwala, S. (2007). A computational framework for simulation of biogeochemical tracers in the ocean. Global Biogeochemical Cycles, 21(3), GB3001. https://doi.org/10.1029/2007GB002923

Khatiwala, S. (2018). Transport Matrix Method software for ocean biogeochemical simulations [Software]. Zenodo. https://doi.org/10.5281/zenodo.1246300

Khatiwala, S. (2023). Fast spin‐up of geochemical tracers in ocean circulation and climate models. Journal of Advances in Modeling Earth Systems, 15(2), e2022MS003447. https://doi.org/10.1029/2022MS003447

Khatiwala, S. (2024). Efficient spin‐up of Earth System Models using sequence acceleration. Science Advances, 10(18), eadn2839. https://doi.org/10.1126/sciadv.adn2839

Khatiwala, S. (2025). tmm4py: Global ocean biogeochemical modeling in Python with the Transport Matrix Method [Data and Software]. https://doi.org/10.5281/zenodo.14878521

Khatiwala, S., Graven, H. D., Payne, S., & Heimbach, P. (2018). Changes to the air‐sea flux and distribution of radiocarbon in the ocean over the 21st century. Geophysical Research Letters, 45(11), 5617–5626. https://doi.org/10.1029/2018GL078172

Khatiwala, S., Schmittner, A., & Muglia, J. (2019). Air‐sea disequilibrium enhances ocean carbon storage during glacial periods. Science Advances, 5(6), eaaw4981. https://doi.org/10.1126/sciadv.aaw4981

Khatiwala, S., Visbeck, M., & Cane, M. (2005). Accelerated simulation of passive tracers in ocean circulation models. Ocean Modelling, 9(1), 51–69. https://doi.org/10.1016/j.ocemod.2004.04.002

Kriest, I., & Oschlies, A. (2015). Mops‐1.0: Towards a model for the regulation of the global oceanic nitrogen budget by marine biogeochemical processes. Geoscientific Model Development, 8(9), 2929–2957. https://doi.org/10.5194/gmd‐8‐2929‐2015

Lam, S. K., Pitrou, A., & Seibert, S. (2015). Numba: A LLVM‐based Python JIT compiler. In Proceedings of the 2nd Workshop on the LLVM Compiler Infrastructure in HPC. Association for Computing Machinery. https://doi.org/10.1145/2833157.2833162

Liang, J.‐H., Deutsch, C., McWilliams, J. C., Baschek, B., Sullivan, P. P., & Chiba, D. (2013). Parameterizing bubble‐mediated air‐sea gas exchange and its effect on ocean ventilation. Global Biogeochemical Cycles, 27(3), 894–905. https://doi.org/10.1002/gbc.20080

McManus, J. F., François, R., Gherardl, J. M., Kelgwin, L., & Drown‐Leger, S. (2004). Collapse and rapid resumption of Atlantic meridional circulation linked to deglacial climate changes. Nature, 428(6985), 834–837. https://doi.org/10.1038/nature02494

Mills, R. T., Adams, M. F., Balay, S., Brown, J., Dener, A., Knepley, M., et al. (2021). Toward performance‐portable PETSc for GPU‐based exascale systems. Parallel Computing, 108, 102831. https://doi.org/10.1016/j.parco.2021.102831

Missiaen, L., Bouttes, N., Roche, D. M., Dutay, J., Quiquet, A., Waelbroeck, C., et al. (2020). Carbon isotopes and Pa/Th response to forced circulation changes: A model perspective. Climate of the Past, 16(3), 867–883. https://doi.org/10.5194/cp‐16‐867‐2020

Nicholson, D. P., Emerson, S. R., Khatiwala, S., & Hamme, R. C. (2011). An inverse approach to estimate bubble‐mediated air‐sea gas flux from inert gas measurements. In S. Komori, W. McGillis, & R. Kurose (Eds.), The 6th International Symposium on Gas Transfer at Water Surfaces (pp. 223–237). Kyoto University Press.

Nicholson, D. P., Khatiwala, S., & Heimbach, P. (2016). Noble gas tracers of ventilation during deep water formation in the Weddell Sea. In The 7th International Symposium on Gas Transfer at Water Surfaces, IOP Conference Series: Earth and Environmental Science. https://doi.org/10.1088/1755‐1315/35/1/012019

Okuta, R., Unno, Y., Nishino, D., Hido, S., & Loomis, C. (2017). CuPy: A NumPy‐Compatible Library for NVIDIA GPU calculations. In Proceedings of Workshop on Machine Learning Systems (LearningSys) in the 31st Annual Conference on Neural Information Processing Systems (NIPS).

Orr, J. C., Dutay, J., Najjar, R. G., Bullister, J., & Brockmann, P. (1999). CFC‐HOWTO: Internal OCMIP report (Technical Report). LSCE/CEA Saclay.

Orr, J. C., Najjar, R., Sabine, C. L., & Joos, F. (1999). Abiotic‐HOWTO. In Internal OCMIP report (p. 15). LSCE/CEA Saclay. Retrieved from www.ipsl.jussieu.fr/OCMIP

Primeau, F. (2005). Characterizing transport between the surface mixed layer and the ocean interior with a forward and adjoint global ocean transport model. Journal of Physical Oceanography, 35(4), 545–564. https://doi.org/10.1175/jpo2699.1

Rempfer, J., Stocker, T. F., Joos, F., Lippold, J., & Jaccard, S. L. (2017). New insights into cycling of 231Pa and 230Th in the Atlantic Ocean. Earth and Planetary Science Letters, 468, 27–37. https://doi.org/10.1016/j.epsl.2017.03.027

Sasaki, Y., Kobayashi, H., & Oka, A. (2022). Global simulation of dissolved 231Pa and 230Th in the ocean and the sedimentary 231Pa/230Th ratios with the ocean general circulation model COCO ver4.0. Geoscientific Model Development, 15(5), 2013–2033. https://doi.org/10.5194/gmd‐15‐2013‐2022

Schmittner, A., & Somes, C. J. (2016). Complementary constraints from carbon (13C) and nitrogen (15N) isotopes on the glacial ocean’s soft‐tissue biological pump. Paleoceanography, 31(6), 669–693. https://doi.org/10.1002/2015PA002905

Séférian, R., Berthet, S., Yool, A., Palmiéri, J., Bopp, L., Tagliabue, A., et al. (2020). Tracking Improvement in simulated marine biogeochemistry between CMIP5 and CMIP6. Current Climate Change Reports, 6(8), 95–119. https://doi.org/10.1007/s40641‐020‐00160‐0

Seltzer, A., Davidson, P. W., Shackleton, S. A., Nicholson, D. P., & Khatiwala, S. (2024). Global ocean cooling of 2.3°C during the Last Glacial Maximum. Geophysical Research Letters, 51(9), e2024GL108866. https://doi.org/10.1029/2024GL108866

Seltzer, A., Nicholson, D. P., Smethie, W. M., Tyne, R. L., Roy, E. L., Stanley, R., et al. (2023). Dissolved gases track deep ocean ventilation processes in the deep North Atlantic. Proceedings of the National Academy of Sciences, 120(11), e2217946120. https://doi.org/10.1073/pnas.2217946120

Siddall, M., Henderson, G. M., Edwards, N. R., Frank, M., Mu, S. A., Stocker, T. F., & Joos, F. (2005). 231Pa/230Th fractionation by ocean transport, biogenic particle flux and particle type. Earth and Planetary Science Letters, 237(1–2), 135–155. https://doi.org/10.1016/j.epsl.2005.05.031

Siddall, M., Khatiwala, S., van de Flierdt, T., Jones, K., Goldstein, S., Hemming, S., & Anderson, R. F. (2008). Towards explaining the Nd paradox using reversible scavenging in an ocean general circulation model. Earth and Planetary Science Letters, 274(3–4), 448–461. https://doi.org/10.1016/j.epsl.2008.07.044

Thiele, G., & Sarmiento, J. L. (1990). Tracer dating and ocean ventilation. Journal of Geophysical Research, 95(C6), 9377–9391. https://doi.org/10.1029/jc095ic06p09377

Vance, D., Little, S. H., de Souza, G. F., Khatiwala, S., Lohan, M. C., & Middag, R. (2017). Silicon and zinc biogeochemical cycles coupled through the Southern Ocean. Nature Geoscience, 10(3), 202–206. https://doi.org/10.1038/ngeo2890

van Hulten, M., Dutay, J., & Roy‐Barman, M. (2018). A global scavenging and circulation ocean model of thorium‐230 and protactinium‐231 with improved particle dynamics (NEMO‐ProThorP 0.1). Geoscientific Model Development, 11(9), 3537–3556. https://doi.org/10.5194/gmd‐11‐3537‐2018

Verdy, A., & Mazloff, M. R. (2017). A data assimilating model for estimating southern ocean biogeochemistry. Journal of Geophysical Research: Oceans, 122(9), 6968–6988. https://doi.org/10.1002/2016JC012650

Wilson, J. D., Andrews, O., Katavouta, A., de Melo Virissimo, F., Death, R. M., Adloff, M., et al. (2022). The biological carbon pump in CMIP6 models: 21st century trends and uncertainties. Proceedings of the National Academy of Sciences, 119(29), e2204369119. https://doi.org/10.1073/pnas.2204369119

Yool, A., Popova, E. E., & Anderson, T. R. (2013). MEDUSA‐2.0: An intermediate complexity biogeochemical model of the marine carbon cycle for climate change and ocean acidification studies. Geoscientific Model Development, 6(5), 1767–1811. https://doi.org/10.5194/gmd‐6‐1767‐2013

Yu, E., Francois, R., & Bacon, M. P. (1996). Similar rates of modern and last‐glacial ocean thermohaline circulation inferred from radiochemical data. Nature, 379(6567), 689–694. https://doi.org/10.1038/379689a0