MPDATA

This article has multiple issues. Please help improve it or discuss these issues on the talk page. (Learn how and when to remove these messages) This article may require cleanup to meet Wikipedia's quality standards. The specific problem is: Reorganising the lead section into a logical order so that the first para really tells us what it is about would help a great deal. Please help improve this article if you can. (August 2026) (Learn how and when to remove this message) This article's lead section may be too long. Please read the length guidelines and help move details into the article's body. (August 2026) (Learn how and when to remove this message) (Learn how and when to remove this message)

In numerical analysis of partial differential equations, MPDATA is a family of iterative finite-difference/finite-volume methods for numerical integration of hyperbolic differential equations modelling conservation laws of the form:

$${\displaystyle \partial _{t}(G\psi )+\nabla \cdot ({\mathbf {v}}\psi )=GR,}$$

1

where ${\displaystyle \psi ({\mathbf {x}},t)}$ is an advected scalar field (or advectee); ${\displaystyle {\mathbf {v}}=G{\dot {\mathbf {x}}}}$ is the flow velocity vector field (or advector), ${\displaystyle G({\mathbf {x}})}$ may either play the role of the fluid density, the Jacobian of coordinate transformation from Cartesian to curvilinear framework ${\displaystyle {\mathbf {x}}}$, or their product; ${\displaystyle R({\mathbf {x}})}$ combines all source terms.

MPDATA stands for Multidimensional Positive Definite Advection Transport Algorithm. The algorithm was formulated by Piotr K. Smolarkiewicz (and is also referred to as Smolarkiewicz's method) at the US National Center for Atmospheric Research - NCAR (at the time, Smolarkiewicz was a fellow of the NCAR Advanced Study Program and a recent graduate from Krzysztof Haman's group at the University of Warsaw; the seminal 1983 and 1984 MPDATA publications mention both institutions).

The crux of the method lies in iterative application of the upwind scheme. The first iteration employs the advective velocity ${\displaystyle {\mathbf {v}}}$, while each subsequent iteration employs a so-called antidiffusive velocity which corrects solution from prior iteration reducing the numerical diffusion. The antidiffusive velocities can be derived through modified equation analysis of the upwind scheme and feature cross-dimensional dependencies (i.e., applying MPDATA in multiple dimensions is not equivalent to application of one-dimensional MPDATA in all dimensions), the scheme is thus not dimensionally split, hence "M" in the algorithm name. Since each iteration of the scheme constitutes a forward-in-time upwind pass, the scheme inherits characteristics of the upwind scheme: CFL stability criterion, conservativeness, embarrassingly parallel domain decomposition, and sign-preservation. For non-negative fields ${\displaystyle \psi }$, sign-preservation translates to positive definiteness, hence the "PD" in the algorithm name. Application of the corrective iterations improves scheme convergence rate compared with first-order upwind. Depending on the MPDATA variant, convergence of up to third-order in time and space can be achieved.

While the original formulation of MPDATA employed structured grids, the algorithm has been subsequently also formulated for unstructured grids. Despite being formulated for and named in reference to advection problems, as any other advection numerical scheme, MPDATA also applies to solutions of advection-diffusion as well as diffusion-only problems if Fickian diffusive terms are expressed as advective fluxes (approach referred to as the pseudo-velocity technique).

Description of the basic scheme in 1D

MPDATA is inherently multi-dimensional, and primarily used in computational fluid dynamics where the advective volocities and problem geometries are variable in time. Still, the key idea underlying the MPDATA approach can be conveyed with a basic example of solenoidal stationary flow in one dimension (i.e., ${\displaystyle {\mathbf {v}}=[u]}$ constant in time and space), without coordinate transformation (${\displaystyle G=1}$), for the case of homogeneous advection (${\displaystyle R=0}$) of a nonnegative scalar field (${\displaystyle \psi \geq 0}$), with the following flux form of the advection equation:

$${\displaystyle \partial _{t}\psi +\partial _{x}(u\psi )=0.}$$

2

Upwind discretisation of the problem on a regular staggered grid with a time step ${\displaystyle \Delta t}$ and a grid step ${\displaystyle \Delta x}$, with ${\displaystyle n=t/\Delta t}$, ${\displaystyle i=x/\Delta x}$, and the half-integer spatial indices corresponding to grid-cell walls:

	${\displaystyle u_{i-3/2}}$		${\displaystyle u_{i-1/2}}$		${\displaystyle u_{i+1/2}}$		${\displaystyle u_{i+3/2}}$
...	|	·	|	·	|	·	|	...
		${\displaystyle \psi _{i-1}}$		${\displaystyle \psi _{i}}$		${\displaystyle \psi _{i+1}}$

can be formulated with:

$${\displaystyle {\frac {\psi _{i}^{n+1}-\psi _{i}^{n}}{\Delta t}}+{\frac {\overbrace {f(\psi _{i}^{n},\psi _{i+1}^{n},u_{i+1/2}^{n})} ^{\text{right-hand wall flux}}-\overbrace {f(\psi _{i-1}^{n},\psi _{i}^{n},u_{i-1/2}^{n})} ^{\text{left-hand wall flux}}}{\Delta x}}=0}$$

3

with the flux function defined using positive and negative parts of ${\displaystyle u_{i\pm 1/2}}$ as:

$${\displaystyle f(\psi _{l},\psi _{r},u)=\!\!\overbrace {\frac {u+|u|}{2}} ^{\text{positive part}}\!\psi _{l}+\!\!\overbrace {\frac {u-|u|}{2}} ^{\text{negative part}}\!\psi _{r}.}$$

4

Introducing the non-dimensional Courant number ${\displaystyle C=u{\Delta t}/{\Delta x}}$, the resultant explicit-in-time scheme (referred to as "upwind", "upstream" or "donor-cell"), for a constant ${\displaystyle C}$ reads:

$${\displaystyle \psi _{i}^{n+1}=\psi _{i}^{n}-\left[{\frac {C+|C|}{2}}(\psi _{i}^{n}-\psi _{i-1}^{n})+{\frac {C-|C|}{2}}(\psi _{i+1}^{n}-\psi _{i}^{n})\right]}$$

5

which is conservative, sing-preserving and stable for ${\displaystyle |C|\leq 1}$. However, the scheme is only first-order accurate in space (${\displaystyle \Delta x}$) and time (${\displaystyle \Delta t}$), and incurs numerical diffusion which can be quantified with the modified-equation analysis by substituting the discretised values other than ${\displaystyle \psi _{i}^{n}}$ in (5) with their Taylor expansions (using the Big O notation):

$${\displaystyle {\begin{aligned}\psi _{i}^{n+1}&=\psi _{i}^{n}+{\frac {\left.\partial _{t}\psi \right|_{i}^{n}}{1!}}(\,+\Delta t\,)+{\frac {\left.\partial _{t}^{2}\psi \right|_{i}^{n}}{2!}}(\,+\Delta t\,)^{2}+O(\Delta t^{3})\\\psi _{i+1}^{n}&=\psi _{i}^{n}+{\frac {\left.\partial _{x}\psi \right|_{i}^{n}}{1!}}(+\Delta x)+{\frac {\left.\partial _{x}^{2}\psi \right|_{i}^{n}}{2!}}(+\Delta x)^{2}+O(\Delta x^{3})\\\psi _{i-1}^{n}&=\psi _{i}^{n}+{\frac {\left.\partial _{x}\psi \right|_{i}^{n}}{1!}}(-\Delta x)+{\frac {\left.\partial _{x}^{2}\psi \right|_{i}^{n}}{2!}}(-\Delta x)^{2}+O(\Delta x^{3})\end{aligned}}}$$

6

which yields (up to second-order terms):

$${\displaystyle {\begin{aligned}\left.\partial _{t}\psi \right|_{i}^{n}\Delta t+\underbrace {\partial _{t}^{2}\psi } _{u^{2}\partial _{x}^{2}\psi }|_{i}^{n}{\frac {\Delta t^{2}}{2}}&={\frac {C}{2}}\left(\psi _{i-1}^{n}-\psi _{i+1}^{n}\right)+{\frac {|C|}{2}}\left(\psi _{i-1}^{n}-2\psi _{i}^{n}+\psi _{i+1}^{n}\right)\\&=-\,C\,\Delta x\,\partial _{x}\left.\psi \right|_{i}^{n}+{\frac {|C|}{2}}\,\Delta x^{2}\,\partial _{x}^{2}\left.\psi \right|_{i}^{n}\end{aligned}}}$$

7

in which the highlighted second-order time derivative can be replaced with second-order spatial derivative using the Cauchy-Kovalevskaya procedure (i.e., by substituting eq. (2) into a time-derivative of eq. (2) giving ${\displaystyle \partial _{t}^{2}\psi =-u\partial _{x}\partial _{t}\psi =u^{2}\partial _{x}^{2}\psi }$) yielding the following, so-called, modified equation:

$${\displaystyle \partial _{t}\psi |_{i}^{n}+u\partial _{x}\psi |_{i}^{n}=\underbrace {\left(|u|{\frac {\Delta x}{2}}-u^{2}{\frac {\Delta t}{2}}\right)} _{k}\,\partial _{x}^{2}\psi |_{i}^{n}}$$

8

which confirms that, as ${\displaystyle \Delta t\rightarrow 0}$ and ${\displaystyle \Delta x\rightarrow 0}$, the upwind scheme approximates eq. (2), but the truncation error has a leading-order Fickian source term with a coefficient of numerical diffusion given by ${\displaystyle k}$.

The concept of MPDATA lies in reversing the effect of numerical diffusion by using the pseudo-velocity technique that allows to express diffusion terms in transport problems using an equation of the same form as (2). To this end, Smolarkiewicz introduced the anti-diffusive pseudo velocity (note the opposite sign stemming from, de facto, integrating backwards in time to reverse the effect of diffusion):

$${\displaystyle \partial _{t}\psi +\partial _{x}{\big (}\!\!\!\!\!\!\!\!\!\!\!\!\underbrace {k{\frac {\partial _{x}\psi }{\psi }}} _{\text{anti-diffusive velocity}}\!\!\!\!\!\!\!\!\!\!\!\!\psi {\big )}=0}$$

9

and proposed discretisation allowing to perform a corrective step using upwind integration with an anti-diffusive Courant number field:

$${\displaystyle C_{i+1/2}^{m+1}={\frac {\Delta t}{\Delta x}}k_{i+1/2}^{m}\left.{\frac {\partial _{x}\psi }{\psi }}\right|_{i+1/2}^{m}\approx {\begin{cases}0&{\text{ if }}\psi _{i+1}^{m}+\psi _{i}^{m}=0\\\left[|C_{i+1/2}^{m}|-(C_{i+1/2}^{m})^{2}\right]{\frac {\psi _{i+1}^{m}-\psi _{i}^{m}}{\psi _{i+1}^{m}+\psi _{i}^{m}}}&{\text{ otherwise }}\end{cases}}}$$

10

where ${\displaystyle m>0}$ numbers algorithm iterations. The first pass of the scheme is an ordinary upwind integration using ${\displaystyle C^{m=1}=C}$ and yielding ${\displaystyle \psi ^{m=1}}$. In the second pass, the values of the anti-diffusive Courant number ${\displaystyle C^{m=2}}$ are computed based on ${\displaystyle C^{m=1}}$ and ${\displaystyle \psi ^{m=1}}$, and used to perform an upwind anti-diffusive pass yielding ${\displaystyle \psi ^{m=2}}$. Note that even for constant-in-${\displaystyle x}$ physical ${\displaystyle C^{m=1}}$, the anti-diffusive ${\displaystyle C^{m=1}}$ varies in ${\displaystyle x}$ (and thus corresponds to a divergent flow). Subsequent iterations are correcting the integration of the anti-diffusive corrections from previous passes. If ${\displaystyle M}$ iterations are used, ${\displaystyle \psi ^{m=M}}$ is used as the result of integration of the advective term over one ${\displaystyle \Delta t}$.

Note that despite a constant flow field, the anti-diffusive velocity varies in space and time and is not solenoidal.

The values of ${\displaystyle \psi _{i-1}^{n}}$ at the left edge of the domain, and ${\displaystyle \psi _{i+1}^{n}}$ at the right edge of the domain are obtained by specifying boundary conditions.

Algorithm variants and techniques used in concert with MPDATA

The basic formulation of the scheme has been extended yielding a multitude of MPDATA variants, including (in chronological order of literature references):

• Divergent-flow correction applicable for non-solenoidal flows (including the antidiffusive velocities).
• Double-Pass Donor Cell (DPDC) in which an analytical approximation of the sum of all corrective iterations is used to define a single corrective pass
• Non-oscillatory variant featuring a flux limiter to ensure monotone solutions
• Infinite-gauge variant for transport of variable-sign fields (e.g., momentum)
• third-order accurate anti-diffusive velocities for a wide range of flow types.

Note that some of the non-basic variants of the algorithm have larger stencils) extent than the original formulation.

Some of the pertinent numerical techniques used together with MPDATA include:

• application to problems featuring pressure-gradient source terms (e.g., Boussinesq approximation (buoyancy) of the Navier-Stokes equations)
• immersed boundary method
• adaptive mesh refinement controlled by antidiffusive velocities.

Applications

Equation (1), which is numerically approximated with MPDATA, is used in modelling a wide range of transport phenomena across different scales and . MPDATA has documented applications in the following domains:

• numerical weather prediction (NWP) systems:
  • compressible dynamical core for the COSMO system
  • nonhydrostatic finite-volume dynamical core for the ECMWF Integrated Forecast System (IFS)
• ocean modelling:
  • passive tracer advection in terrain-following coordinates in the Regional Ocean Modeling System (ROMS)
• large eddy simulations:
  • moist atmospheric flows (including representation of cloud and precipitation processes)
  • multi-phase particle-laden flows (e.g., contrail formation, aerosol-cloud-precipitation interactions)
  • turbulent magnetohydrodynamics (including Hall magnetohydrodynamics)
  • simulations of wind-driven sand dune dynamics
  • simulations of flows past complex structures (e.g. the Pentagon building and the Palace of Culture and Science in Warsaw)
• direct numerical simulations:
  • multi-phase flows (fog formation)
• solvers for shallow water equations resolving:
  • double-gyre circulation
  • wetting–drying areas
• aerosol size-spectrum dynamics:
  • advection problem in particle size dimension
  • two-dimensional advection in size-spectral and spatial dimensions
• mathematical finance:
  • solutions to the Black-Scholes-type equations (a stochastic differential equation transformed into a partial differential equation using Itô's lemma) cast as an advective transport problem for pricing of European and American options

Open-source implementations

• reusable MPDATA libraries and packages:
  • libmpdata++: structured-grid MPDATA implemented in C++ using Blitz++ for array handling and array-valued expressions, and OpenMP for multi-threading. API and usage examples documented in a journal article. Featuring distributed-memory parallelism using Message Passing Interface. Supports integration in 1-, 2- and 3-dimensional domains. Implements arbitrary number of corrective iterations, the infinite-gauge, non-oscillatory and fully-third-order variants. Features implementations of a shallow-water equation system solver, and a conjugate-residual Poisson equation solver for handling pressure terms in momentum conservation equations. The library has been developed as a dynamical core for the UWLCM large eddy simulation system. Released under the GNU GPL v3 license.
  • PyMPDATA structured-grid MPDATA implemented in Python using NumPy and Numba for array handling, just-in-time compilation and multi-threading, supporting distributed-memory parallelism via Numba-MPI. Supports integration in 1-, 2- and 3-dimensional domains, arbitrary number of corrective iterations, implements infinite-gauge, non-oscillatory and double-pass-donor-cell algorithm variants. Pure-Python wheels available on PyPI, release tarballs archived at Zenodo. Released under the GNU GPL v3 license together with a set of examples in a form of Jupyter notebooks. Web-based code-coupled documentation in English.
• MPDATA implementations integrated in other software:
  • PISM (Parallel Ice Sheet Model, GPL v3 license) includes implementation of the non-oscillatory variant of two-dimensional MPDATA in C++
  • AtmosFOAM (a derivative of OpenFOAM; OpenFOAM is distributed under the GPL v3 license) features implementations of MPDATA in both C++ and Python/Numba
  • mpdat_2d: a FORTRAN 77 subroutine implementing 2-dimensional structured-grid MPDATA including the non-oscillatory option (likely a derivative of the original MPDATA implementation created at NCAR) has been publicly released as a part of babyEULAG code. No documentation or README, unspecified license and authorship.
  • smolar: structured-grid MPDATA implemented as a C library supporting multi-threading and multi-process parallelism, as well as offloading of anti-diffusive velocity computation to GPU using CUDA. Code publicly available on GitHub without release versioning, unspecified license, documentation in Russian language credited to Russian Academy of Sciences.
  • advectionHPCtester: structured-grid 3-dimensional implementation of the basic scheme in FORTRAN. Code archived on Zenodo. Released under the GNU LGPL.

References