The non-conservation of potential vorticity by a dynamical core compared with the effects of parametrized physical processes

Numerical models of the atmosphere combine a dynamical core, which approximates solutions to the adiabatic, frictionless governing equations for ﬂuid dynamics, with tendencies arising from the parametrization of other physical processes. Since potential vorticity (PV) is conserved following ﬂuid ﬂow in adiabatic, frictionless circumstances, it is possible to isolate the effects of non-conservative processes by accumulating PV changes in an air-mass-relative framework. This ‘PV tracer technique’ is used to accumulate separately the effects on PV of each of the different non-conservative processes represented in a numerical model of the atmosphere. Dynamical cores are not exactly conservative because they introduce, explicitly or implicitly, some level of dissipation and adjustment of prognostic model variables which acts to modify PV. Here, the PV tracers technique is extendedtodiagnosethecumulativeeffectofthenon-conservationofPVbyadynamicalcore anditscharacteristicsrelativetothePVmodiﬁcationbyparametrizedphysicalprocesses. QuantiﬁcationusingtheMetOfﬁceUniﬁedModelrevealsthatthemagnitudeofthenon-conservationofPVbythedynamicalcoreiscomparabletothosefromphysicalprocesses. Moreover,theresidualofthePVbudget,whentracingtheeffectsofthedynamicalcore andphysicalprocesses,isatleastanorderofmagnitudesmallerthanthePVtracers associatedwiththemostactivephysicalprocesses.Theimplicationofthisworkisthat thenon-conservationofPVbyadynamicalcorecanbeassessedincase-studieswitha fullsuiteofphysicsparametrizationsanddirectlycomparedwiththePVmodiﬁcationby parametrizedphysicalprocesses.Thenon-conservationofPVbythedynamicalcoreis showntomovethepositionoftheextratropicaltropopausewhiletheparametrizedphysical processeshavealessereffectatthetropopauselevel.


Introduction
Potential Vorticity (PV) thinking has become a key concept in dynamical meteorology.PV has two key properties, conservation (Ertel 1942) and invertibility as discussed in Hoskins et al. (1985).Conservation means that, in the absence of diabatic and frictional processes, PV is advected like a tracer.Invertibility means that the PV distribution, with appropriate boundary conditions, is sufficient to diagnose all of the dry dynamical variables to the approximation of a given balance condition.The usefulness of PV thinking depends on the accuracy of the balanced dynamics.Davis et al. (1996) demonstrated that most of the dynamics of an intense extratropical cyclone could be quantified using the balance equations of Charney (1955).McIntyre and Norton (2000) showed that higher order PV-based balanced models were capable of producing simulations "remarkably similar" to the full unbalanced equations for shallow-water simulations.
Considering the PV conservation in a numerical model of the atmosphere Davis et al. (1993) partitioned PV into a set of tracer diagnostics to explicitly integrate the cumulative effects of parametrized physical processes in a study of cyclogenesis.Combined with the piecewise PV inversion method of Davis and Emanuel (1991) this allowed them to assess the impact of non-conservative processes on a cyclone's circulation.The PV diagnostics did have limitations.Davis et al. (1993) noted differences between the PV tracers and the PV diagnosed from model variables and attributed this to numerical truncation errors in updating PV.Stoelinga (1996) discussed this difference in more detail and attributed it to the difference between the explicit PV integration and the model dynamics which are not designed to conserve PV exactly.Zhang et al. (2008) demonstrated that inconsistencies between tracer advection and the dynamical core of an atmospheric model can produce significant biases in modelling chemical transport.Whitehead et al. (2015) assessed the consistency of several dynamical cores with their respective tracer advection schemes by using potential vorticity in an idealised baroclinic wave test.In this test there is no diabatic heating or friction, so the only source of PV non-conservation is dissipation, implicit or explicit, induced by the dynamical core.Whitehead et al. (2015) tested the consistency of this dissipation with the dissipation induced by the tracer advection scheme for a range of dynamical cores and demonstrated that each dynamical core produces values of PV inconsistent with the tracer advection scheme but with structure and amplitude differing between dynamical cores.Whitehead et al. (2015) used the consistency between a dynamical PV and a tracer of PV to rank the different dynamical cores.
Tracers of PV have been used for two things: assessing the dynamical impacts of parametrized physical processes and diagnosing inconsistencies between dynamical cores and tracer advection schemes.In this study we demonstrate that, by extending the PV diagnostics of Davis et al. (1993), we can explicitly diagnose the inconsistency in PV (defined by Whitehead et al. (2015)) in a simulation which also has parametrized physical processes.We show that this "dynamics-tracer inconsistency" is comparable to the effects on PV of parametrized physical processes for a case study with the Met Office's Unified Model (MetUM) and that the majority of the "dynamicstracer inconsistency" in the MetUM can be attributed to nonconservation of PV by the dynamical core.
We will consider the effects of the nonconservation of PV by the dynamical core in the context of previous studies with the PV tracers in the MetUM.Using the PV tracers, Chagnon et al. (2013) identified a dipole in diabatically modified PV across the extratropical tropopause for a case study of an extratropical cyclone.The extratropical tropopause corresponds to a sharp gradient in PV and can be defined as a surface of PV with 2 PVU (PV units) (Hoskins et al. 1985;Reed 1955).The tropopause dipole consisted of positive diabatically generated PV in the stratosphere and negative diabatically generated PV in the troposphere either side of the 2 PVU surface, suggesting a sharpening of the tropopause PV gradient due to non-conservative processes (Chagnon et al. 2013).Chagnon and Gray (2015) showed that a tropopause dipole was also present in three other case studies of extratropical cyclones.We will show that the nonconservation of PV by the dynamical core of the MetUM acts to diminish the tropopause dipole in our case study.
The structure of the paper is as follows.Section 2 describes the methods of analysis as well as introducing the case study.Section 3 presents the results obtained from a detailed investigation of the PV budget in the MetUM and the implication for the tropopause dipole.Section 4 provides a summary and discussion of the results.A large amount of specific notation is used in this paper so a notation table is provided in appendix A.

Methodology
In this section we describe the methods used in this study.A simulation was performed using the MetUM version 7.3.Details of the MetUM are provided below with a description of the PV tracers method and the case study used.Offline trajectory calculations have also been used in this study and details of the method are provided at the end of this section.

Met Office Unified Model
The MetUM is an operational numerical weather prediction model (Davies et al. 2005).The dynamical core of the MetUM approximates a two time level, semi-implicit, semi-Lagrangian solution to the governing equations of the atmosphere.The prognostic variables in the MetUM are wind (u = (u, v, w)), potential temperature (θ), the mixing ratios of moisture variables (specific humidity, cloud ice and cloud liquid), density (ρ) and Exner pressure (Π).The variables in the MetUM are placed on an Arakawa C-grid with Charney-Phillips staggering in the vertical.
The MetUM contains various parametrizations to account for physical processes that are either not resolved or not represented within the dynamical core: radiation (Edwards and Slingo 1996); microphysics (Wilson and Ballard 1999); non-orographic gravitywave drag (Scaife et al. 2002) and orographic gravity-wave drag (Webster et al. 2003); convection (Gregory and Rowntree 1990); and turbulent mixing represented by the boundary-layer scheme of Lock et al. (2000).
The simulation in this paper uses a previously operational limited area domain known as the NAE (North Atlantic and Europe) domain nested within a global domain using the method described in Davies (2013); Fig. 2 shows the extent of the NAE domain.The lateral boundary conditions used were produced from the operational runs of the global model and the start dump used is from the operational NAE analysis.The NAE has 0.11°horizontal grid spacing (approximately 12 km) with 70 stretched vertical model levels up to 80 km.The model levels in the MetUM use a terrain following hybrid-height coordinate that gradually flattens at higher altitudes (Davies et al. 2005).We use the standard timestep for this domain of 5 minutes.
A semi-Lagrangian method does not explicitly conserve any variables.However, the MetUM incorporates an Eulerian discretisation of the continuity equation to ensure local mass conservation (Davies et al. 2005).No explicit diffusion is applied to prognostic variables; the diffusion in the MetUM is entirely implicit and a result of the cubic interpolation (quintic for moisture variables) used in the semi-Lagrangian scheme.The MetUM contains a modified vertical interpolation for θ, described in section 6 of Davies et al. (2005), that is required for stability.This modified vertical interpolation is applied up to a height of 3.4 km for our simulation.
The MetUM solves the governing equations of the atmosphere using a "predictor-corrector" method; an initial "predictor" is made of the prognostic variables at time level n + 1 and is refined using a set of "correctors".The full method and governing equations are set out by Davies et al. (2005) and the inclusion of parametrized physical processes to the equations is presented in Diamantakis et al. (2007).We present a simplified description of the method, based on section 5 of Davies et al. (2005) and section 2 of Diamantakis et al. (2007), so that the budget of PV in the MetUM can be described precisely in the following section.
The discretisation of the governing equations which the MetUM aims to approximate is given by where X is a vector of the prognostic variables in the MetUM; L and N are the linear and nonlinear dynamics terms respectively; SP and FP are the tendencies of the "slow" and "fast" parametrized physical processes; a subscript "d" denotes evaluation at departure points in the MetUM's semi-Lagrangian method; α is a time-weighting coefficient (typically 0.7); n and n + 1 are the time levels; and ∆t is the time step.Figure 1 shows a schematic of a single timestep of the MetUM which demonstrates the application of the predictor-corrector method to solving Eq. ( 1).The components of this figure are now described.First, a set of increments due to the slow physical processes (microphysics (mic), radiation (rad) and gravity-wave drag (gwd)) are calculated from the set of prognostic variables at the start of the timestep.The increment to prognostic variables due to slow physical processes can be written as (2) This article is protected by copyright.All rights reserved.Next solutions to the thermodynamic, moisture and momentum equations are approximated.This is the predictor step for X.Equation ( 1) is solved explicitly with time level n + 1 values replaced with time level n estimates.This can be written as (3) where X (1) is the first estimate of X n+1 .Note that the increments due to slow physical processes, interpolated to departure points, are added on at this stage.This predictor step could be split into two stages with an intermediate state of where X sl is the set of prognostic variables after the semi-Lagrangian dynamics only and the slow physical processes increments act as the first corrector for X.
The next corrector adds the effects of fast physical processes (convection (con) and boundary layer (bl)) using the most up to date estimates of the prognostic variables.The fast physical processes are calculated sequentially for stability.This can be written as (5) At this stage density and Exner pressure have not been updated from their time-level n estimates.The continuity equation is discretised in an Eulerian form and Exner pressure is used to couple the prognostic variables using the ideal gas equation of state.The back substitution of the equations to replace time-level n estimates with n + 1 values leads to a Helmholtz-type equation to solve; the full equations are given in appendix B of Davies et al. (2005).This is known as the pressure solver and is written as another corrector: where N * is the latest estimate of N.
At the end of the timestep the MetUM modifies the prognostic variables in clouds to eliminate supersaturation and account for the additional latent heat release (known as cloud balancing).This can be considered as a final physical process corrector:

PV Tracers
A set of PV tracers are integrated online in the MetUM.The method is similar to that developed by Davis et al. (1993) and was first applied to the MetUM by Gray (2006).The method works by partitioning and integrating PV.The evolution of PV, in Lagrangian form is given by (Ertel 1942), where q = 1 ρ (∇ × u + 2Ω) • ∇θ is PV; θ and F are the diabatic heating and friction respectively; and Ω is Earth's rotation rate vector.
The general method is to integrate Eq. ( 8) along trajectories over a forecast of time T : where S is the right-hand side of Eq. ( 8) and is partitioned into different physical processes (S = S i ) resulting in a set of PV tracers (q i ) from the integration of S i starting with each q i = 0.
The PV diagnostics are essentially mimicking the behaviour of the numerical weather prediction model in terms of PV, allowing the tendencies of each parametrized physical process to be partitioned and accumulated separately.There is an implicit assumption here that the effects of each parametrized physical process can be separated.In practice the PV tracers will often have large cancelling terms between compensating processes.It is important to consider all terms in the PV budget to assess where this is the case.
We partition PV into an advection-only PV tracer (q adv ) and a set of physics PV tracers ( q phys ).Each PV tracer, apart from the advection-only PV tracer, is set to zero everywhere at the initial time and the advection-only PV tracer is initialised as equal to the PV diagnosed from the prognostic variables X (diagnosed PV).At the lateral boundaries of a limited area domain each PV tracer, apart from the advection-only PV tracer, is set to zero and the advection-only PV is set equal to the diagnosed PV, at each timestep.This is because there is no prior information on the history of the air parcels at the lateral boundaries so they are treated like initial conditions.
In the MetUM tracers are advected by the flow resolved in the model using its semi-Lagrangian advection scheme.Passive tracers can also be transported by the effects of sub-grid parametrizations for turbulence and convection.The turbulence scheme acts primarily within the planetary boundary layer and acts as a down-gradient mixing where the diffusion coefficient is dependent upon the Richardson number (Lock et al. 2000).In the simulations shown here the turbulent diffusion acts in the vertical only.The convection scheme is based on a mass flux formulation (Gregory and Rowntree 1990) and passive tracer transport depends upon the diagnosed profile of convective updraught and downdraught mass fluxes.The parametrized subgrid scale motions have no horizontal component across the sides of a grid-box and, by construction, the updraughts and downdraughts are exactly compensated by vertical motion in the remainder of the column (above the grid box on the Earth's surface) such that the area-averaged vertical motion at each level equals the resolved vertical motion in that grid box.
Although PV is materially conserved in adiabatic, frictionless flows, Haynes and McIntyre (1987) showed that it behaves very This article is protected by copyright.All rights reserved.

Accepted Article
differently to a passive tracer once non-conservative processes act.They showed that there can be no net transport of PV across any isentropic surface, even when diabatic and frictional processes are acting.They called this remarkable property the "PV impermeability theorem" (Haynes and McIntyre 1990) and it can be regarded as a consequence of Kelvin's circulation theorem.It applies generally to fully compressible non-hydrostatic dynamics and implies that the cross-isentropic flux of PV must be identically zero.This applies to the resolved motions in a model and notional sub-grid transports.Furthermore, the theorem also holds for other vertical coordinates.For example, in pressure coordinates there can be no net transport of the vertical component of vorticity across pressure levels.Therefore, since the sub-grid turbulence and convective schemes describe only vertical fluxes and the nontransport of vorticity must hold, these transport schemes are not applied to the PV tracers and only advection by the resolved 3-D motion transports them.The tracer advection scheme in the MetUM also has the option to apply a conservation correction and a monotone correction.However, neither is applied to the PV tracers.The tracer advection scheme simply amounts to updating a variable with its departure point value, obtained by interpolating the tracer at time-level n to the departure point of the trajectory calculated by the semi-Lagrangian scheme of the MetUM.In this way we are calculating the left hand side of Eq. ( 8), partitioned by each PV tracer, at each timestep, following the resolved flow of the forecast.
Although there are no PV fluxes across isentropic surfaces, diabatic and frictional effects have an important influence on the PV distribution via diabatic mass fluxes across isentropic surfaces, horizontal divergent flow and frictional torques.Each PV tracer, apart from the advection-only PV tracer, accumulates increments in PV due a specific parametrized physical process at each timestep.The PV increment is calculated as the difference in PV before and after adding the increments to the prognostic variables.In this way we are calculating the right hand side of Eq. ( 8), partitioned by each parametrized physical process, at each timestep.
The slow physical processes are calculated in parallel and each increment is calculated independently using X n .The increments in PV due to slow physical processes are calculated as ∆qsp = q(X n + ∆Xsp) − q(X n ), (11) where q(•) means a calculation of PV as a function of the given variables in the argument and sp is microphysics, radiation or gravity-wave drag.The fast physical processes are calculated sequentially so the PV increments are calculated as ∆qcon = q(X (1) + ∆Xcon) − q(X (1) ), (12) for convection and ∆q bl = q(X (1) + ∆Xcon + ∆X bl ) − q(X (1) + ∆Xcon), (13) for the boundary layer scheme, where q(X (2) ) ≡ q(X (1) + ∆Xcon + ∆X bl ).The increment from cloud balancing is also included with the physics PV tracers, calculated as ∆q cloud = q(X (4) ) − q(X (3) ).

(14)
Therefore the equations for updating the PV tracers are given by for the advection-only PV, where the "d" subscript indicates that the tracer is evaluated at departure points in the MetUM's semi-Lagrangian scheme; for slow physical processes; and for fast physical processes, where f p is convection, boundary layer or cloud balancing.
In the previous equations where q(•) is evaluated we use a modified version of the MetUM's standard diagnostic PV calculation.PV is calculated at the corners of grid points, on model levels where ρ is stored, such that the calculation of the vertical component of vorticity requires no averaging.Each component of PV is calculated using centred differences with prognostic variables averaged when required.The result is then linearly interpolated to the centres of the grid points in the horizontal then linearly interpolated to model levels where θ is stored in the vertical.The modification is that we have also included vertical velocity components in the calculation of PV which improved the PV budget presented later in this paper.

Case Study
A 36 h simulation was performed using PV tracers with the MetUM version 7.3 initialised from Met Office operational analysis at 1200 UTC 28th November 2011.This period corresponds to a case study from the DIAMET (DIAbatic influences on Mesoscale structures in ExTratropical storms) project (Vaughan et al. 2015) of a double cold front passing over the UK.Diabatic heating and cooling rates in the cold front due to phase changes of water were investigated by Dearden et al. (2014) using in-situ flight data for this case study.More in depth details of the meteorological conditions during the case study are given in Dearden et al. (2014).
Our case study was also included in Chagnon and Gray (2015) (case II) as another example of the tropopause dipole of Chagnon et al. (2013) and corresponds to the case with the weakest tropopause dipole, specifically less diabatically-generated positive PV in the stratosphere.We have investigated this case study to explain the lack of diabatically-generated positive PV in the stratosphere using the PV tracer technique.We investigate the effects of the non-conservation of PV by the dynamical core on the tropopause dipole.

Trajectories
PV tracers tell us about the behaviour of physical processes following the resolved flow of the forecast.However tracer advection schemes on a grid inevitably involve dissipation as features cascade to scales too small to represent.Trajectories calculated from model winds also tell us about processes following the resolved flow of the forecast but do not include any dissipation.We use the trajectory calculation method of Wernli and Davies (1997) (LAGRANTO, Sprenger and Wernli (2015)) with hourly wind data output from the model to provide an alternative view of the forecast flow.
Trajectories are used to calculate PV redistributed by advection as an alternative to the advection-only PV tracer.The advectiononly PV tracer at any time tells us the PV that the resolved air parcel had at the start of the model run.PV calculated from trajectories also tells us the PV that the resolved air parcel had at the start of the model run but is not affected by the numerical diffusion of the tracer advection scheme at each timestep.However, it is dependent on the consistency of the trajectory calculations with the MetUM.The trajectory calculations will differ from the model trajectories because of the different numerical schemes used and the lower frequency of data used for the offline calculation of trajectories (hourly).
This article is protected by copyright.All rights reserved.

Results
If all processes modifying PV were accounted for by the PV tracers then the accumulated effects of parametrized physical processes ( q phys ) would be equivalent to the total change in PV given by the diagnosed PV minus the advection-only PV (q − q adv ).We define ε as the difference between these two measures of PV change, such that where q n ≡ q(X n ) is the diagnosed PV.Note that ε is the same as qr in Stoelinga (1996).
Figure 2b shows ε and Fig. 2a shows the sum of physics PV tracers for comparison, each at the end of the 36 h forecast interpolated to the 320 K isentropic surface.The field of ε has structure and amplitude that is comparable to the sum of physics PV tracers and therefore represents an important contribution to the PV budget.

PV Budget
In this section, by considering the evolution of the PV budget across a timestep, we show that ε is completely accounted for by three terms: "dynamics-tracer inconsistency" (∆ε I ) to be defined below based on the inconsistency investigated by Whitehead et al. (2015); "missing PV" (∆ε M ) which accounts for any increments in PV not attributed to a dynamical or physical process; and a "splitting error" (∆ε S ) which accounts for the difference between numerical diffusion acting on multiple tracers of PV and the numerical diffusion acting on a single field representing the sum of those PV tracers.In this section we define these three terms mathematically and their numerical form in the MetUM.In principle a PV budget for any numerical model of the atmosphere could be closed using just these three terms.
To calculate a closed PV budget we must account for all changes in PV across a single timestep, such that the PV at time level n + 1 is equal to the PV at time level n plus all the changes in PV in that time-step.Considering every change in PV across a single timestep of the MetUM (as in Fig. 1) we can write: We can attribute each of the terms in Eq. ( 19) to increments in PV related to dynamical and physical processes.The first increment in PV in Eq. ( 19) (q(X sl ) − q(X n )) is a result of the semi-Lagrangian dynamics, which we can define as ∆q sl , such that q(X sl ) − q(X n ) = ∆q sl . (20) The next increment in PV in Eq. ( 19) (q(X (1) ) − q(X sl )) is a result of adding the increments to prognostic variables from the slow physical processes to the latest estimate of X.This is approximately equivalent to updating the slow physics PV tracers: where qsp is the set of PV tracers for slow physical processes.The equation is not exact because of the nonlinearity associated with the calculation of the PV increments due to slow physical processes (Eq.( 11)) in parallel and the order in which the increments are added.The next increment in PV in Eq. ( 19), is the increment in PV due to the fast physical processes (Eqs ( 12) and ( 13)).The next increment in PV in Eq. ( 19) (q(X (3) ) − q(X (2) )) is a result of the pressure solver, which we can define as ∆q solver , such that q(X (3) ) − q(X (2) ) = ∆q solver . (23) The final increment in PV in Eq. ( 19), is the increment in PV due to cloud balancing (Eq.( 14)).We can define the increments due to fast physical processes as The increment due to cloud balancing ∆q cloud has been grouped with fast physical processes for convenience.Whitehead et al. (2015) define the inconsistency between a dynamical core and tracer advection as the difference between the evolution of PV calculating by integrating the governing equations in a dynamical core and the advection of a tracer of PV.We can define this "dynamics-tracer inconsistency" for a single timestep as ∆ε I by comparing PV obtained by solving the adiabatic and frictionless governing equations (q n + ∆q sl ) with PV advected using the tracer advection scheme (q n d ), such that Ideally the increment in PV due to the pressure solver from Eq. ( 23) would be included because the pressure solver involves the solution of the continuity equation and the backsubstitution to complete the solution of the thermodynamic and momentum equations.However, the pressure solver also couples the parametrized physics to the dynamics so it is not completely attributable to adiabatic and frictionless dynamics.With all increments in PV described in terms of dynamics and physics, we can use Eq. 26 to rewrite Eq. ( 19) as where ∆ε M is the "missing PV" and includes the increment in PV due to the pressure solver as well as accounting for the nonlinearity in the calculations of PV increments due to slow physical processes from Eq. ( 21).Rearranging Eq. ( 18) for time level n + 1 gives and we can eliminate the terms describing PV tracers in Eq. ( 27) using the definitions of PV tracer updates (Eqs (15), ( 16) and ( 17)).However we first need to account for the difference between numerical diffusion acting on multiple tracers of PV and the numerical diffusion acting on a single field representing the sum of those PV tracers highlighted by the placement of the d subscript in Eq. ( 29).This "splitting error" is a result of diffusion in the tracer advection scheme which is entirely implicit for the operational MetUM and in the simulation performed here.We define the "splitting error" as This article is protected by copyright.All rights reserved.
Accepted Article We can now eliminate all terms describing PV tracers in Eq. ( 27) using Eqs ( 28), ( 29) and ( 30) and the definitions given by Eqs ( 15), ( 16) and ( 17), which gives the result Since ε is by definition zero everywhere at the start of a forecast, Eq. ( 31) tells us that the gap in the PV budget can only be due to the accumulation of the three terms defined in this section: "dynamics-tracer inconsistency" (∆ε I ), "missing PV" (∆ε M ) and a "splitting error" (∆ε S ) and modifications due to the implicit diffusion of ε from carrying it forward on departure points in Eq. ( 31).

Dynamics-Tracer Inconsistency
The previous derivation allows us to describe ε completely as the accumulation of three terms: dynamics-tracer inconsistency; missing changes in PV across a single timestep; and the difference between advecting a single tracer of PV and advecting multiple tracers of PV.In this section we will show that dynamics-tracer inconsistency is the dominant contribution to ε.
The dynamics-tracer inconsistency is calculated in the MetUM at each timestep using Eq. ( 26).It is then accumulated in the same way as a PV tracer: The residual of the PV budget is calculated as and can only be due to the "missing PV" and the "splitting error" (i.e. the time integral of ∆ε M + ∆ε S ). Figure 2 shows the integrated dynamics-tracer inconsistency (ε I ) and the residual PV (εr) for a 36 hour forecast.Dynamics-tracer inconsistency accounts for most of ε and the residual PV is generally more than an order of magnitude smaller than ε and has smallerscale structure.Therefore the pressure solver and nonlinearities in calculations of PV increments (∆ε M ) and the "splitting error" (∆ε S ) are comparatively small contributions to the PV budget for our case study.

Non-conservation of PV by the dynamical core
In this section we show that the dominant process contributing to dynamics-tracer inconsistency is the non-conservation of PV by the dynamical core rather than the numerical dissipation of a PV tracer in the tracer advection scheme.We have already shown that ∆ε S is a small contribution to the PV budget which tells us that the numerical diffusion acting on multiple tracers of PV is approximately the same as the numerical diffusion acting on a single field representing the sum of those PV tracers.In this section we show that the numerical dissipation of a single tracer of PV is small compared to the dynamics-tracer inconsistency.What is the true change in PV integrated along the resolved flow of the forecast?Dynamics-tracer inconsistency arises both from non-conservation of PV by the dynamical core and the numerical dissipation in the tracer advection scheme.If the tracer advection scheme were perfectly conservative, the total change in PV following an air-mass would be the difference between the diagnosed PV and the PV from the origin of the trajectory that the air-mass followed through the whole forecast (q − q origin ).
This article is protected by copyright.All rights reserved.

Accepted Article
Including the PV at the origin of trajectories in the PV budget (Eq.( 18)) gives which highlights a notional partition between non-conservation of PV by the dynamical core (the left bracket) and non-conservation associated with numerical dissipation in the tracer advection scheme acting on the advection-only PV tracer (right bracket).
To estimate the relative magnitudes of these two contributions, trajectories were released from a regular grid on the 500 hPa surface at the end of the forecast (T+36) and calculated backwards in time to the start of the forecast using hourly 3-D wind output from the MetUM (Section 2.4).The initial PV field is interpolated to each trajectory location at t = 0, providing an estimate of q origin which can be associated with the grid point that the trajectory "arrives on" at the end of the forecast.This technique is called a "reverse domain filling" (RDF) trajectory calculation since a map of q origin (x) is obtained (Fig. 3c).
Figure 3 compares three different measures of PV: the diagnosed PV (q), the advection-only PV (q adv ), and the PV calculated from RDF trajectories (q origin ).The fields are shown at 500 hPa because the back trajectory calculations were initialised on pressure levels.The 500 hPa surface is found using linear interpolation while assuming a logarithmic variation of pressure with height consistent with other MetUM diagnostics.Owing to the exact conservation implied by the RDF technique, the maxima and minima in the field are given by the extrema in the initial PV distribution (at the origin locations which in general are not at 500hPa due to vertical motion).Fine scales are generated through stirring by advection and there is no dissipation in the reverse domain filling calculation to smooth the smallscale structure.In contrast the advection-only tracer experiences numerical diffusion.This acts to remove the smallest structures and to fill in some regions with intermediate PV values (for example, around the cyclonic spiral to the southwest of Iceland in Fig. 3).The highest PV values in RDF PV over southern England are also reduced in the PV tracer, presumably by mixing in the tracer calculation.In contrast, the diagnosed PV (from the prognostic variables) shows much lower values than q adv or q origin within the low PV air to the south of Iceland.Part of this difference is associated with physical processes and part with the non-conservation by the dynamical core.
Figure 4 shows ε calculated from Eq. ( 18) and ε calculated from the first bracket in Eq. (34) using the RDF estimate q origin .There are considerable differences between the two terms, mainly in terms of fine-scale structure fluctuating about zero.As already discussed the fine scale structure arises from lack of dissipation in the RDF calculation and also small errors associated with the offline calculation of long trajectories used in the RDF calculation.However, it can be seen that ε calculated from RDF trajectories accounts for most of the larger scale and magnitude PV anomalies seen in ε calculated from Eq. ( 18).Therefore, we can conclude that numerical diffusion of the advection-only PV tracer is not the major contribution to the dynamics-tracer inconsistency in ε.
Tracer advection within the MetUM has the option to run with various different interpolation schemes.Running the same simulation while varying the interpolation scheme used for the PV tracers from linear to quintic makes very little difference to the dynamics-tracer inconsistency (not shown).This suggests that the numerical diffusion of the PV tracers associated with the interpolation to departure points in the semi-Lagrangian advection scheme is not the major contribution to the dynamics-tracer inconsistency.
Having eliminated other options, the conclusion is that the tracer advection used for the PV tracers is more conservative, in terms of PV, than the dynamical core.The majority of the dynamics-tracer inconsistency must therefore be due to the nonconservation of PV by the dynamical core.
How should we interpret this non-conservation of PV by the dynamical core?Like PV, θ is conserved in the absence of diabatic and frictional processes.Therefore it can be partitioned into a set of tracers in the same way as PV (Martínez-Alvarado and Plant 2014) and will also have an associated dynamics-tracer inconsistency.However, θ is a prognostic variable in the MetUM so the changes to θ in the dynamical core are essentially identical to tracer advection (depending on the interpolation schemes used).The dynamics-tracer inconsistency for θ (calculated with Eq. ( 26)) would therefore be close to zero and tests have shown it is close to zero (Martínez-Alvarado, personal communication).
If we were to run the same study as Whitehead et al. (2015) with the MetUM we would get different answers depending on whether we use PV or θ.The prognostic variable θ tells us that the MetUM has a consistent dynamical core and tracer advection scheme.However, the diagnostic variable PV shows a large dynamics-tracer inconsistency.PV is diagnosed as a function of the gradients of the prognostic variables, each of which are updated separately (Eq.( 1)).Therefore the Lagrangian equation for PV (Eq.( 8)) is not respected exactly by the dynamical core of the MetUM.The dissipation of prognostic variables θ and u in the dynamical core will act like the effects of heating and friction terms on PV.
We could also associate the non-conservation of PV by the dynamical core to an implicit representation of physical processes at small scales.Kunkel et al. (2014) used a PV tracer in an adiabatic and frictionless simulation to assess the impact of inertia-gravity waves near the tropopause.Kunkel et al. (2014) found systematic differences between the diagnosed PV and tracer PV in the regions of inertia-gravity waves.Small-scale physical processes would act to modify PV but should not be expected to modify passive tracers in the same way.

Tropopause Dipole
Using the PV tracers method Chagnon et al. (2013) showed that the accumulated effects of parametrized physical processes contributed to a sharpening of the tropopause PV gradient for a case study of an extratropical cyclone.In this section we investigate the effects of the non-conservation of PV by the dynamical core on the tropopause PV gradient for our case study.Chagnon et al. (2013) showed, for their case study, that the 2 PVU surface of the advection-only PV tracer coincided with the 2 PVU surface of the diagnosed PV.This meant that the direct modifications of PV by non-conservative processes did not act to change the position of the tropopause.They showed that by summing the values of the near tropopause total change in PV (q − q adv ), binned by the advection-only PV tracer, that the average change in PV for values of advection-only PV less than 2 PVU (tropospheric) was negative and the average change in PV for values of advection-only PV greater than 2 PVU (stratospheric) was positive with a zero value at 2 PVU (Fig. 6 in Chagnon et al. (2013)).Chagnon and Gray (2015) repeated this diagnostic for three more extratropical cyclones.They showed that our case study (case II in Chagnon and Gray (2015)) contains regions of tropopause sharpening but the average strength of the tropopause dipole is weaker and the dipole does not have a strong positive stratospheric PV modification.
Figure 5a shows the same diagnostic as  Figure 4.The difference between the accumulated effects of parametrized physical processes ( q phys ) and the total change in PV calculated using (a) advection-only PV (q adv ) and (b) PV traced back along a trajectory (qorigin), at 500 hPa and T+35 h so that we do not miss any shallow PV anomalies that may be important.We also include only gridpoints within 2.5 km in the vertical of the tropopause, excluding the boundary layer using the MetUM's diagnosis of boundary layer height.The mean mass of the included gridpoints is 3.18 × 10 10 Kg, with a maximum of 4.21 × 10 10 Kg and a minimum of 5.69 × 10 9 Kg.Whitehead et al. (2015) found large amounts of dynamicstracer inconsistency where isentropes intersected the ground.We find the same result for our case study.However, this low level "dynamics-tracer" inconsistency is found to largely cancel out tendencies between the radiation and boundary layer parametrization schemes such that the total change in PV (q − q adv ) is much smaller in magnitude.By excluding the boundary layer we avoid seeing large signals of opposite sign in Fig. 5a.
Figure 5a does show a dipole in the total change in PV.The dipole is consistent with the one in Chagnon and Gray (2015) with the addition of a positive PV anomaly at q adv between 6 and 8 PVU which can be attributed to the use of many vertical levels.More relevant to this study is the large difference between the sum of physics PV tracers and the total change in PV (q − q adv ).This difference is mostly accounted for by the non-conservation of PV by the dynamical core (ε I ).In Fig. 2 there is a difference between the position of the 2 PVU contour of the advection-only PV and the 2 PVU contour of the diagnosed PV that is directly caused by This article is protected by copyright.All rights reserved.Figure 6 shows the same diagnostic as Fig. 5 as a function of lead time for the accumulated effects of parametrized physical processes and the total change in PV.The difference between Fig. 6a and Fig. 6b can be attributed to the non-conservation of PV by the dynamical core ( I ).The sum of physics PV tracers shows a clear dipole at a value of advection only PV tracer close to the 2 PVU tropopause; however this is not present for the total change in PV which appears to have a faint dipole that drifts with time.This suggests that the non-conservation of PV by the dynamical core acts to move the 2 PVU tropopause position whereas the parametrized physical processes do not.

Conclusions
A new diagnostic framework is introduced to calculate the nonconservation of PV by the dynamical core of a numerical model of the atmosphere when simulating a realistic case study with a full suite of physics parametrizations.The non-conservation of PV by the dynamical core has been considered in the context of PV tracers based on the method introduced by Davis et al. (1993).Whitehead et al. (2015) used tracers of PV to diagnose inconsistencies between dynamical cores and tracer advection schemes but applied to idealised simulations without any parametrization of physical processes.A "dynamics-tracer inconsistency" diagnostic has been incorporated into the PV tracers method in the MetUM and used to diagnose the nonconservation of PV by the dynamical core.We have shown that, for our case study, the non-conservation of PV by the dynamical core has a comparable contribution to the PV budget to that of parametrized physical processes.Similar results have also been produced with the dynamics-tracer inconsistency diagnostic in other case studies (not shown).
Discrepancies between the PV diagnosed from the prognostic variables of a model and the PV tracers have been previously noted.Davis et al. (1993) attributed the difference to numerical errors in the explicit integration of PV.Stoelinga (1996) attributed the difference to using a numerical model that does not conserve PV explicitly.Gray (2006) and Chagnon et al. (2013) attributed the difference to the amplified effects of diffusion across multiple tracers.In reality all of these terms could be important and will have differing importance for different numerical models.We have introduced a framework that can account for each of these effects separately and we can calculate the relative importance of these terms for the PV tracers method applied to any numerical weather prediction model.
The residual in the PV budget is generally more than an order of magnitude smaller than the dominant physical processes when the non-conservation of PV by the dynamical core is accounted for.Currently the largest part of the residual in the PV budget comes from the pressure solver.If the residual in the PV budget were larger then a method to sensibly partition the PV increment from the pressure solver would need to be developed for the PV tracers.A possible method to include the pressure solver would be to run two timesteps of the model in parallel: one regular timestep and one adiabatic and frictionless timestep.The latter would be used to calculate the dynamics PV increments in calculating the dynamics-tracer inconsistency.This method is generic for any dynamical core and would also be a closer match to the inconsistency defined by Whitehead et al. (2015).It has not been attempted for this study because the dynamics-tracer inconsistency diagnosed from the semi-Lagrangian dynamics step was the major term missing from the PV budget.
The version of the MetUM (7.3) used in this study does not add any explicit diffusion to the PV tracers.In the initial formulation of the PV tracers method Davis et al. (1993) chose to add diffusion to the PV increments at each timestep to mimic the effect of explicit thermal diffusion on small-scale anomalies.Our approach is to use the PV tracers to assess the behaviour of the numerical model itself.
An advantage of the PV tracers is that the cumulative effects of different processes can be quantified in a single simulation, as opposed to model sensitivity studies where the different experiments will in general have different model trajectories.This is important because the processes interact nonlinearly and so each process depends sensitively upon the model trajectory.The technique could be used to relate forecast errors to the processes contributing and to identify systematic model error.
It has been shown that numerical weather prediction models systematically smooth the tropopause PV gradient with lead time (Gray et al. 2014).Gray et al. (2014) hypothesised that the smoothing of the PV gradient was due to an underrepresentation of diabatic processes consistent with the development of a diabatic PV dipole shown by Chagnon et al. (2013).However we have shown that the the non-conservation of PV by the dynamical core has a strong effect on the tropopause and could also explain the smoothing of the PV gradient.The results of Chagnon et al. (2013) and Chagnon and Gray (2015) do implicitly include the non-conservation of PV by the dynamical core because they look at differences between the diagnosed and advection-only PV.However by attributing PV to dynamics-tracer inconsistency This article is protected by copyright.All rights reserved.the accumulated effects of parametrized physical processes ( q phys ) and (b) the total change in PV (q − q adv ).The solid and dashed lines show positive and negative anomalies respectively.

Accepted Article
we have reduced the uncertainty in the PV budget and therefore reduced the uncertainty in the individual PV tracers, further validating the approach in Chagnon et al. (2013) of looking at the effects of individual physical processes on the PV dipole with the caveat that the non-conservation of PV by the dynamical core should also be considered.By looking at the evolution of numerical solutions to idealised cases of frontogenesis past the point of frontal collapse, Visram et al. (2014) suggested that insufficient Lagrangian conservation of PV can cause a degradation to the long-term solutions of forecasts.The nonconservation of PV by the dynamical core would be a direct cause of this.However, we refrain from describing the nonconservation of PV by the dynamical core as "model error" because it is necessary to have some form of dissipation in numerical models of the atmosphere and this dissipation may also be linked to unrepresented small-scale physical processes.By diagnosing non-conservation of PV by the dynamical core with the dynamics-tracer inconsistency diagnostic we can assess the Lagrangian conservation of PV in case studies and differentiate between physical processes and model error.

Figure 1 .
Figure 1.A schematic of a single timestep of the MetUM.

Figure 2 .
Figure 2. (a) The sum of the physics PV tracers ( q phys ).(b) ε as defined by Eq. (18).(c) Dynamics-tracer inconsistency.(d) The residual in the PV budget (εr in Eq. (33)).All plots are linearly interpolated to the 320 K isentrope for a 36 h forecast.The thick black line in each plot is the 2 PVU line of the diagnosed PV and the dashed line is the 2 PVU line of the advection-only PV tracer Fig. 6 in Chagnon et al. (2013) but integrated over many vertical levels by weighting the gridpoints by mass rather than an area average over individual vertical levels, with Fig 5b showing the total mass associated with each bin.The diagnostic was integrated over many vertical levels This article is protected by copyright.All rights reserved.

Figure 5 .
Figure 5. (a) The near-tropopause mass-weighted average PV in bins of 0.25 PVU of the advection-only PV tracer.(b) The total mass in each 0.25 PVU bin of advectiononly PV.Results shown for an integration time of 36 hours.

Figure 6 .
Figure6.The near-tropopause mass-weighted average PV in bins of 0.25 PVU of the advection-only PV tracer against time (hours since forecast initialisation) for (a) the accumulated effects of parametrized physical processes ( q phys ) and (b) the total change in PV (q − q adv ).The solid and dashed lines show positive and negative anomalies respectively.