Global mass fixer algorithms for conservative tracer transport in the ECMWF model

Various mass fixer algorithms (MFAs) have been implemented in the Integrated Forecasting System (IFS) of the European Centre for Medium-Range Weather Forecasts (ECMWF) to ensure mass conservation of atmospheric tracers within the semi-Lagrangian (SL) advection scheme. Emphasis has been placed in implementing schemes that despite being primarily global in nature adjust the solution mostly in regions where the advected field has large gradients and therefore interpolation (transport) error is assumed larger. The MFAs have been tested in weather forecast, idealised and atmospheric dispersion cases. Applying these fixers to specific humidity and cloud fields did not change the accuracy of 10-day forecasts. In other words, global mass tracer conservation is achieved without deteriorating the solution accuracy. However, for longer forecast timescales or for forecasts in which correlated species are transported, experiments suggest that MFAs may improve IFS forecasts.


Introduction
A drawback of semi-Lagrangian (SL) transport schemes, such as the one used by the European Centre for Medium-Range Weather Forecasts (ECMWF) Integrated Forecasting System (IFS Ritchie et al., 1995) is that they do not formally conserve mass as the pointwise nature of the SL method does not take into account grid-box size and fluxes.Between the beginning and the end of each time step, the total model mass can differ by a very small amount.This difference, although not significant for the timescales of numerical weather prediction (NWP), may accumulate in the long run.A systematic drift in the total mass of air (or a tracer field) will eventually affect the quality of the forecast (Thuburn, 2008).
As NWP models become more complex, the number of tracers increases and therefore the requirement for conservative schemes becomes more important.Furthermore, as the resolution increases towards cloud-resolving scales it becomes increasingly desirable from the parametrisation point of view to have a mass-conserving advection scheme as this may improve further the simulation of cloud processes.
SL advection (SLA) consists of two steps which do not -in principle -ensure conservation of mass: (i) finding departure points and (ii) interpolating the advected field to the departure point location.However, the choice of method for (i) and (ii) has a considerable impact for the amount of the mass non-conservation.
There is a class of SL schemes, the so-called inherently conserving schemes, which are able to achieve global, local and consistent mass conservation for tracer and air-mass fields.Two examples are the SLICE (semi-Lagrangian inherently conserving and efficient) transport scheme (see Zerroukat and Allen, 2012) and CSLAM (conservative semi-Lagrangian multi-tracer) transport scheme (see Lauritzen et al., 2010).These schemes are an application of a finitevolume-type discretisation approach on the semi-Lagrangian continuity equation.In general, they are complex algorithms difficult to implement efficiently in an existing operational model which uses a "traditional" SL method.Although inherently conserving SL methods are not currently used in weather forecasting operations there are schemes in this family which are competitive or even more efficient than their Eulerian finite-volume conservative counterparts for applications where a large number of tracers is advected (multitracer simulations).CSLAM is an example of such a method while another example of a recent development based on the Published by Copernicus Publications on behalf of the European Geosciences Union.
An alternative low-computational-cost approach to ensure global mass conservation which can be easily applied on traditional SL methods is the mass fixer algorithm (MFA).The task of a MFA is to change the tracer concentrations after SLA in such way that the mass before and after advection is the same.A general problem of MFAs is to identify regions where it is most appropriate to change the solution of the SL scheme.
Different MFAs implement different strategies for distributing the global mass loss or gain.The simplest ones correct the solution uniformly by simply scaling each grid-point value with the ratio of the global mass before and after advection.This approach is currently used in IFS when long time integrations take place in order to correct the total model mass and that of long-lived tracers (Flemming and Huijnen, 2011).
More sophisticated MFAs attempt to compute a correction which is proportional to the smoothness of the solution.A larger correction is applied in areas where the solution has large gradients and therefore the error is larger, and a very small correction where the solution is smooth and the error is small.
The aim of the paper is to present tracer MFAs that were recently implemented in IFS in model cycle 39r1.Using this model cycle as the base for our experiments we shall discuss results from NWP forecasts, long-range forecasts where the mass fixers are applied to humidity and cloud fields as well as idealised tracer and volcanic plume forecasts.Availability of globally mass-conserving schemes for tracers can be an important addition to IFS-based prediction systems such as the EC-Earth (Hazeleger et al., 2012) climate model or atmospheric composition forecast systems where aerosols, greenhouse and reactive gases are transported (Hollingsworth et al., 2008).
The paper is structured as follows.The amount of the nonconservation by the SL advection scheme of IFS is demonstrated in Sect. 2. Section 3 describes the implemented MFA.Their impact on the simulated fields in different applications is discussed in Sect. 4. Conclusions are presented in Sect. 5.

Air and tracer global mass conservation in IFS
In a 10-day IFS forecast, at the current operational resolution T1279L137 (approximately 16 km in grid-point space on 137 levels) using a 10 min time step, the total model air mass increases by less than 0.01 % of its initial value.The formulation of the continuity equation, based on the Ritchie and Tanguay (1996) scheme (see also ECMWF, 2012, Sect. 3.6.2),plays an important role in achieving this accuracy.Orography is removed from the advected mass field resulting in a much smoother field which can be accurately interpolated to the Lagrangian grid (departure points).
Global conservation errors in tracer advection are larger and depend on the smoothness of the field.For example, smoother fields such as ozone and specific humidity have smaller conservation errors than fields with sharp features such as cloud fields.This is demonstrated in Fig. 1 where the global mass conservation error is displayed for ozone, specific humidity (Q), liquid cloud water content (CLWC), and cloud ice water content (CIWC) for the same number of time steps (1440) at different resolutions using two approaches for the interpolation to the departure point.Mass conservation is represented in Fig. 1 by a line identical to the horizontal 0 axis.The global mass conservation error for a tracer φ is expressed as a percentage of its initial mass: , where M φ 0 and M φ t are the initial-and current-step global tracer mass.
In the forecast experiments of Fig. 1 all parametrisations of sink and source terms have been switched off.This allows one to test the performance of the advection scheme using real orography.In addition, the following two interpolation methods have been used: (i) the quasi-cubic ECMWF interpolation (Ritchie et al., 1995) with a quasi-monotone limiter and (ii) a linear interpolation (indicated with LIN in plots).Method (i) is used in IFS operationally for Q and ozone while method (ii) is used operationally for the rougher cloud fields.The experiments are run at the following horizontal and vertical resolutions: (i) T159 L60 i.e.T159 in the horizontal (approximately equal to 125 km) with 60 levels in the vertical, (ii) T159 L91, (iii) T1279 L91 (approximately 16 km in the horizontal), and (iv) T1279 L137.To allow direct comparisons of the mass conservation error per time step, the four forecasts in Fig. 1 have been run for the same number of time steps.At coarse horizontal resolution (T159) the time step is six times longer (60 min) than the corresponding time step for high resolution (T1279).
The results shown in Fig. 1 indicate that the global mass conservation error per time step tends to decrease as the resolution increases.However, when horizontal resolution is increased from T159 to T1279, the accumulated error at t = 10 days decreases only for CLWC, CIWC with cubic interpolation while remains roughly the same for the remaining fields.It seems that the opposite is true when vertical resolution increases, the accumulated error at t = 10 days decreases except for CLWC and CIWC with cubic interpolation.So there are differences between interpolation schemes and between fields of different smoothness but the overall indication is that in the IFS system mass conservation of tracers tends to improve globally as resolution increases and the best way to demonstrate this is by comparing Fig. 1a with d.

Description of the MFAs
The transport problem we consider here is the advection of a scalar field φ χ which represents the mass mixing ratio of a tracer: where ρ χ and ρ are the tracer and air density respectively and S represents sources or sinks that may be present.Consider SL time stepping from t to t + t: , where d denotes the departure point computed by the trajectory algorithm and φ t χ ,d is obtained by interpolating the known field φ t χ at the computed departure point.If S = 0 then the global volume integral of ρφ χ at t and t + t (on the model grid) should not change as this represents the total mass of χ and the only process operating is advection (transport).However, in practice, as the interpolation scheme generates errors this global conservation law is violated.
Global MFAs of different sophistication are described in the published literature for SL transport models.In general, any MFA will compute the global tracer mass immediately before and after the advection step.Then a small correction is computed for each grid point in such a way that this global error is eliminated.In the simplest version of the proportional or multiplicative fixer of Rasch and Williamson (1990), each grid-point value is multiplied by the ratio of the mass before and after advection.Here, we will focus on the more local algorithms.In particular, the following algorithms will be discussed: (i) the quasi-monotone Bermejo and Conde (2002) scheme, (ii) Zerroukat (2010) scheme, (iii) the quasi-monotone Priestley (1993) scheme, and (iv) the Mc-Gregor (2005) scheme.These algorithms have been implemented in IFS and will be summarised in the following paragraphs.It should be noted that their implementation is threedimensional (3-D) given that semi-Lagrangian advection in IFS is fully three-dimensional.
To describe these different fixers, as implemented in IFS, we use the following notation: K is the number of model levels, starting from the top of the atmosphere and ending on the surface.Each model level has N grid points.Each grid box has horizontal surface area A j and height z j k where z j k denotes the height of the j th model grid point of the k th level.The total mass of a tracer χ with mass mixing ratio φ = ρ χ /ρ, where ρ is the air-density field, is The hydrostatic approximation (valid in IFS) p = −ρg z has been used in Eq. ( 2) to eliminate z.
During the advection step, a tracer field φ 0 (i.e. the field before the advection step takes place) is interpolated to the departure point field (Lagrangian grid) and changes to φ * while its total mass changes from M 0 to M * : Use of p * j k in Eq. ( 3) reflects the change of the surface pressure field due to advection.MFAs aim to correct φ * so that a new field is derived which has a total mass equal to M 0 .

Bermejo and Conde (BC) scheme
The Bermejo and Conde (2002) algorithm is derived by a variational principle.It computes a new quasi-monotone field minimising its distance from the original one subject to the constraint of global mass conservation.The correction added at each grid point depends on an estimate of the interpolation error.The global norm of this correction field has the smallest possible magnitude that can give mass conservation and monotonicity.In the original publication, the scheme was tested on idealised two-dimensional cases of advection.Here it has been implemented in IFS in 3-D mode and has been tested on active meteorological fields.
Let φ 1 be the field which minimises the square of the weighted norm: min where w j k is a non-negative weighting factor.Having w j k = 0 means that the corresponding grid-point value is not altered and is not included in the cost function.A solution to Eq. ( 4) is found using a Lagrange multiplier approach.The cost function is defined seeking a pair of values φ 1 , λ such that Solving these two equations we obtain ) where the weight w j k depends on the solution smoothness.
We choose it to be proportional to the difference between the quasi-cubic, quasi-monotone interpolated field φ * and the linear one φ L : The above weights are used to compute a "local correction", i.e. the global mass surplus or deficit is distributed unevenly to different grid points depending on the smoothness of the solution which is measured by the difference between a highand a low-order interpolant.For the IFS implementation, β was set to 1 as tests showed no benefit from using the recommended value β = 3.In fact, higher values led to sharper, bigger size increments which may not be desirable for the model stability.
For convenience, in sections that follow, this scheme will be called BC fixer.

Zerroukat's (ZE) scheme
The BC fixer in IFS can also be run in a mode that corresponds to a version of the Zerroukat (2010) fixer.This leads to smoother correction fields.The drawback is that quasi-monotonicity or positive-definiteness cannot be guaranteed.Here an implementation of this scheme is presented which uses the same measure to assess the solution smoothness as the BC scheme, i.e. the difference between a highorder scheme (cubic Lagrange interpolation) and a loworder scheme (linear interpolation).Here, this scheme will be called ZE fixer.It corrects each grid-point value as follows: where M 0 and M * are defined by Eq. ( 3) and again β = 1 is sufficient for practical purposes.If holds then global mass conservation is guaranteed: It is worth noticing that Eq. ( 7) can be re-written in a form that resembles Eq. ( 5): ) This implies that the derived field φ 1 is also a solution of the minimisation problem of Eq. ( 4).One difference between Eq. ( 5) and Eq. ( 8) is the construction of the weights w j k .
Using the unlimited w j k = |φ * j k − φ L j k | β means that all gridpoints will be corrected.The sign of the increment is determined by the sign of δM (which determines the sign of λ): for δM > 0 (surplus) φ 1 j k ≤ φ j k ∀j, k and for δM < 0 (deficit) φ 1 j k ≥ φ j k ∀j, k.However, as this one-directional correction is not limited as in the BC case, it is possible that a new minimum or maximum value may be generated.In practice, if a quasi-monotone scheme was used for advection this happened in less than 0.5 % for humidity grid points but it can sometimes go up to 5 % of grid points for a non-smooth field.Priestley (1993) produced a well-known mass-fixing scheme.Its objective is to compute a globally conserving monotone solution by blending the original high-order solution with a low-order solution thereby departing as little as possible from the high-order one.This is equivalent to finding the highest possible values for the weights α j k such that the "blended" field

Priestley's (PR) scheme
where φ 0 , j, k denotes the set of φ-field values before advection at grid-points surrounding the (j, k) departure point and φ * and φ L the cubically and linearly interpolated fields at the departure point respectively.The two conditions in Eq. ( 9) ensure conservation and monotonicity.The requirement for "highest possible" α values is an accuracy requirement.It ensures that the final solution is as close as possible to the original high-order interpolation field.In regions where the solution is smooth the blended scheme is weighted towards the higher-order solution while in regions with low degree of smoothness it is blended towards the linear solution.
A more detailed step-by-step algorithmic description of Priestley's algorithm is given in the Appendix of Gravel and Staniforth (1994).Priestley's scheme is an iterative scheme.Two options have been implemented: the standard algorithm which will be called here PR and a variant of it, namely PRqm.The latter is essentially the same algorithm, the only difference here is that a quasi-monotone (qm) filter (Bermejo and Staniforth, 1992) has been applied immediately before the application of the fixer.The result of this modification is that the algorithm converges faster.Regardless which variant is used the solution will be always quasi-monotone, the difference is only in the starting values.

Mc Gregor's (MG) scheme
McGregor (2005) scheme which shall be called here MG fixer, is a MFA used in the climate model C-CAM (Conformal-Cubic Atmospheric Model).This is a model using a SL scheme for horizontal advection and a total variation diminishing (TVD) scheme for the vertical advection.MG fixer can be applied to any interpolation technique including linear as opposed to the fixers considered so far which both require that the field is advected using a high-order interpolant.An additional advantage of this scheme is that it is computationally very cheap.However, it does not guarantee monotonicity but only positive definiteness.Furthermore, it differs from the other algorithms presented here, as it does not use a local smoothness criterion to assess how much to correct at each grid-point.At each time step it computes a global diagnostic which judges the overall ability of the advection scheme to accurately advect fields.Nevertheless it does not correct by the same proportion each grid point but is using instead two different scaling factors: one for points that have positive advective increments and one for points that have negative advective increments.It tends to amplify the solution when there is damping and suppress when there is amplification.
The algorithm can be described as follows: -Step 1: compute total mass before and after advection, M 0 and M * as in Eq. ( 3). - Step 2: let a minimum allowed value φ min .Scan each grid point, compute and store: where .
-Step 3: compute total positive and negative increments and their ratio: -Step 4: set α φ = min r, √ r and update: .
The last step is equivalent to and implies that the increment is scaled by a factor α φ which reduces positive increments when their total mass exceeds the total mass of the negative increments.When the opposite is true then positive increments will be amplified and negatives will reduce in magnitude.The new field satisfies the global mass conservation constraint:

The quasi-monotone limiter
The quasi-monotone limiter renders the interpolation locally monotone, i.e. in the vicinity of the departure point the interpolation curve (or multidimensional surface) passing from the departure point field value and the field values of points surrounding the departure point does not generate new minima or maxima.For the tests presented in the following section, two forms of the quasi-monotone Bermejo and Staniforth (1992) mini-max (minimum-maximum) limiter for cubic interpolation will be used: (i) The "default" limiter or filter used operationally in IFS: the scheme is applied immediately after each 1-D cubic interpolation (in longitude, latitude and height) takes place.So, the steps taken are to interpolate in longitude and then apply a 1-D limiter on the interpolated field.Repeat this action for each of the remaining two interpolations (in latitude and height).For brevity this scheme will be called DEF limiter or filter.
(ii) The standard Bermejo and Staniforth (1992) limiter: this shall be called BS limiter or filter.In this case the limiter is applied after all three interpolations have finished, i.e. this is limiting in 3-D at once.
We should also clarify that the term "cubic interpolation" will imply here the quasi-tricubic interpolation scheme used by IFS (linear interpolation along the edges of the stencil, fully cubic in the interior; see Ritchie et al., 1995).

Testing of MFAs in IFS
In Fig. 2 the global conservation error during the advection of Q and CLWC with and without MFA is displayed.It is shown there that application of a MFA eliminates this error.This forecast run has the operational resolution T1279 horizontal with 137 levels and is identical to the one that corresponds to the results of Fig. 1; i.e. there are no sources or sinks of tracer mass.For brevity we display only results from the BC and PR schemes but the other MFAs also give a globally massconserving solution.The mass-conservation error before and after the advection was always close to machine precision.The impact of the BC MFA on Q is demonstrated in Fig. 3. Cubic interpolation is used for the advection of this field.Here, physical parametrisations have been switched on and the setup is the same as in an operational forecast.A single time step increment from the fixer, at t = 24 h and at a model level which over flat terrain is near the 700 hPa pressure level, is compared with the field itself.The figure shows that the computed increments are at least three orders of magnitude smaller than their corresponding field magnitude.The sign is negative due to the fact that at this stage of the forecast, advection increases mass and the fixer has to remove a global surplus.The fixer is acting mainly on areas where large gradients are present where interpolation is expected to be less accurate.In areas where the field is smooth the correction is very small regardless of the field magnitude.Similar results have been produced from runs with the remaining MFAs.For brevity these will not be displayed here but they are publicly available (see Fig. 5 in Sect. 4 in Diamantakis and Flemming, 2013).A zonally and 24 h time-averaged vertical cross section for Q is compared with corresponding cross sections of increment diagnostics in Fig. 4. The average increment is 4-5 orders of magnitude smaller than the magnitude of the field itself.It is concentrated in areas where large amounts of hu-midity are present.It is interesting to notice how similar the zonally and time-averaged increments are for BC, ZE and PRqm.The fact that their difference is small means that the different algorithms converge roughly to the same solution.Larger differences can be noticed when any of the previous three fixers is compared with MG and even larger with PR.
Usually, increments computed by PR differ in sign and magnitude from the other fixers (see also Figs. 4,5).This is because this algorithm computes a quasi-monotone and conservative solution iteratively starting from a cubic interpolated field.In the tests presented here it usually takes 3-4 iterations for PR to converge.During this iterative process both positive and negative increments will be computed to derive a locally monotone solution.Mass has to be removed from overshooting points (negative increment) and added at undershooting points (positive increment).This is not the case with PRqm which starts with a quasi-monotone field having no undershooting or overshooting points and therefore the only action that the algorithm needs to take is to restore global mass conservation.Regarding the remaining fixers it is worth mentioning that (i) ZE produces the smallest, in magnitude, increments but these are slightly more widespread, (ii) BC and PRqm are similar, and (iii) MG produces slightly different patterns than the previous two fixers.
As expected, the quasi-monotone schemes did not produce any overshoots or undershoots.A very small percentage of undershoots (< 0.01 % of total points) was found with MG but no negative values were created.This percentage was larger in the ZE fixer for the cloud fields, slightly exceeding 1.5 %, while it was of similar magnitude for Q (≈ 0.01 %).
Most of these undershoots generated negative values.
In the plots presented here specific humidity was chosen to examine the local behaviour of MFAs.This choice was made due to the meteorological importance of this tracer field and given that it includes regions that are relatively smooth as well as regions with large gradients.The MFA applied to the rougher cloud fields CLWC and CLIC resulted in similar local patterns as shown for Q.The CLWC increments were used as a diagnostic for demonstrating the step by step behaviour of the MFAs.This is shown in Fig. 5 where the scaled global rms and max norms of the of CLWC fixer (absolute) increments are displayed.These are scaled to be the fraction (percentage) of the rms global norm of the advected CLWC field which is representative to its mean value.The plot shows that the smallest increments are computed by the ZE fixer, followed by BC and PRqm while as expected and explained before PR computes the largest increments.MG increments are in the middle between PR and ZE.
As expected PR is the most expensive and MG the cheapest.All algorithms have been parallelised using MPI and open MP directives.

Impact of humidity MFAs on temperature fields in long runs
As there is a strong interaction between humidity and temperature, typically because of radiative effects and cloud microphysics, we shall test in this section to what extent the mass fixer increments on humidity and cloud fields alter the temperature field.To show the impact we carried out four 12-month forecasts with full physics at T159 L137 resolution.This is a standard test of IFS which is done to evaluate whether a new scheme impacts the model's climate.The experiments run are described in Table 1.In Fig. 6a the temperature bias is plotted, i.e. the difference of the vertical cross section of a zonally averaged annual mean temperature field (averaged across the four forecasts) from its corresponding field from the ERA-Interim (ECMWF Reanalysis) run.This figure displays a common problem in semi-Lagrangian models, the extratropical tropopause/lower stratosphere cold bias (see Stenke et al., 2008).For the remaining plots, the difference of the same field (zonally averaged annual mean temperature) from the control run is used.This is done to clearly demonstrate the impact of the changes.As a general rule, warming around the extratropical tropopause (in the region where the blue area in Fig. 6a   Results show that none of these fixers deteriorates an existing cold bias.When the fixers are combined with the DEF limiter the difference is small (results show a marginal improvement and have not been included here).On the contrary a noticeable improvement, i.e. a reduction of the cold bias, can be noticed when they are combined with the BS limiter.This shows in Fig. 6b-f.Good results are obtained with the quasi-monotone algorithms PR and BC.As condition in Eq. ( 9) shows, the PR fixer is limiting the solution using a scheme similar to BS limiter.Bigger positive impact is obtained by fixers that do not guarantee quasi-monotonicity: ZE followed by MG.However, the former generates negatives especially in the cloud fields which are rougher (3-5 % of grid points become negative after correction is applied).This is not the case for the latter where a negative fixer is built in.

Impact on NWP scores in 10-day forecasts
The accuracy of 10-day forecasts is typically assessed using measures that describe the realism of the global geopotential or temperature fields.The forecast fields are compared against the Analysis of the fields and expressed as root mean square error (RMSE) or anomaly correlation coefficient (ACC) (Wilks, 2011).
In general, it is not expected that a global MFA will improve forecasting skill in the short or medium range but neither it should deteriorate the skill.To investigate this the MFAs have been tested running 37 forecast cases, each starting 10 days apart from 01/12/2011 until 25/11/2012.The resolution used is T511 L137 and each forecast is run for 10 days using operational options for the model dynamics and physics.All fixers were activated on Q, CLWC, CIWC, CRWC, and CSWC.Although these tests are specific on moist physics tracers, they do have a general value.We can indirectly measure the impact a fixer has on advection by measuring the overall forecast skill of the experiment: forecast skill deterioration would imply that the tested algorithm deteriorates the accuracy of the advection scheme and there-fore is deemed not suitable for tracer advection.Neutral scores should indicate that the fixer is making the interpolation conservative without damaging solution accuracy at least on the large scale.
Overall, geopotential, wind, and temperature verification scores in the three global regions (Northern Hemisphere, tropics, and Southern Hemisphere) from runs with MFAs are neutral and there is no forecast that is better in terms of ACC and RMSE.An exception is the temperature RMSE in the tropics at upper tropospheric levels which increases up to 0.07 K (from approximately 1.26 to 1.33 K) at t = 10 days when any MFA is applied for humidity and cloud fields with cubic interpolation options.The fixer contributes further (by a small amount) to the existing cold bias.This happens because a small amount of humidity is removed from the atmosphere as a small humidity surplus is detected by the fixer.Reducing the humidity content of the troposphere has in general a cooling effect while the opposite is true for the stratosphere due to reduction of radiative cooling.However, there is no impact on the corresponding ACC scores which remain neutral.

Simulation of correlated tracers
Mass conservation is an important property for atmospheric applications where chemical species are transported.It is also important that existing functional relationships in their concentration are maintained by the advection scheme (see Lauritzen and Thuburn, 2012).The ability of IFS and the newly developed fixers to preserve such relationships has been tested using case 11 from DCMIP (Dynamical Core Model Intercomparison Project; see Ulrich et al., 2012).This is a three-dimensional, passive advection, deformational flow idealised test case in which four tracers are transported.The initial concentration of the first two tracer fields q 1 and q 2 obeys the non-linear relationship: where λ, θ, z is the longitude, latitude and height of a tracer.The first one (q 1 ) is represented by two cosine bells placed at the same height and latitude but at different longitudes.
Results for this test case from IFS runs at T159 horizontal resolution and 137 levels in the vertical (this is close to the recommended resolution for this problem) are plotted in Fig. 7.These plots are correlation plots for the pair (q 1 , q 2 ) at t = 6 days after the initial time.This is half the time required for the tracers to return to their original position; i.e. complete one full rotation around Earth.The initial concentration of these tracers is given by the parabolic dash-dotted black curve.Pairs (q 1 , q 2 ) (red dots) that fall outside the re-gion marked by the dashed-dotted convex shape correspond to unphysical mixing ratios.Real mixing in the atmosphere can only move scatter points to the concave side of the preexisting functional curve along mixing lines (Lauritzen and Thuburn, 2012).Lack of spread indicates that the scheme is overdiffusive as peak values are damped.
The plots show that semi-Lagrangian transport with linear interpolation is excessively diffusive but does not produce any unphysical mixing.The opposite is true when cubic Lagrange is used.It results in a relatively large amount of unphysical mixing and overshoots/undershoots (new maxima/minima are created corresponding to values above 1 and .q 1 -q 2 (xy axes) scatter plots for correlated tracers at t = 6 days.Scatter points (q 1 , q 2 ) at t = 0 follow the upper (parabolic) black dashed-dotted curve.
elative mass residual in volcanic plume simulations (SO 2 ) for different schemes ).Significant improvements can be noticed when a quasi-monotone limiter is used.The DEF limiter, being more strict (and damping) has bigger impact as all points stay inside the convex shape.However, maximum field values are damped.The BS limiter reduces but does not eliminate completely the unphysical mixing occurring with cubic interpolation.However, it preserves better the maxima.When a MFA is combined with the DEF limiter it does not change the mixing further: it preserves equally well the existing tracer correlations as shown in Fig. 7 (compare panels  c and d).It also results in a small further reduction of maximum field values (result not included here).When the fixers are combined with the BS limiter we obtain very similar results with respect to tracer correlations compared with corresponding results from the DEF limiter but slightly improved results in terms of accuracy (preservation of maxima).In this case BC and PR give the best results.They both preserve reasonably well the initial correlation (better than the corresponding run without fixer) and maximum field values are not too far from the analytical values.ZE and MG fixers are not as effective in preserving the functional relationship (especially the latter) as a small proportion of points are outside Experiments defined as in Table 1.

27
Figure 9.Comparison of volcanic plume simulation with and without a mass fixer using quasi-monotone cubic Lagrange at T1279 L91 resolution.The plotted quantity is the total SO 2 content (in kg m −2 ) per model grid-point column.Experiments defined as in Table 1.
the bounded sector.The former can produce small negative values in some regions.But they are both better in preserving the maxima.
In conclusion, applying any of the MFAs did not deteriorate the mixing properties of the advection scheme and in some occasions improved them (e.g.compare Fig. 7e and f).This is a desirable result and suggests that MFAs can be a beneficial addition for a semi-Lagrangian scheme used for transport of chemical tracers.The combination of a MFA with the BS limiter works better and BC and PR seem to give the best results.

Volcanic plume case study
MFAs have also been tested on volcanic plume advection cases.Here a test case is presented where a tracer (SO 2 ) is emitted into the atmosphere by a single point source and then transported by the winds.This case resembles the Grímsvötn volcanic eruption (see Flemming and Inness, 2013).Due to the highly localised nature of the advected plume, this case is a good test for assessing the local behaviour of a global MFA.The striking fact in this simulation is that the plume's total mass is largely overestimated.A conservation error of almost 20 % of the total mass of the field occurs during the first time steps which eventually results in a more than 50 % gain.This is shown in Fig. 8.The greatly improved performance in terms of conservation of the non-limited cubic Lagrange without MFA shown in this plot is due to the presence of large negative undershoots which offset the overshoots when the global integral is computed and is therefore misleading.
Applying a MFA results in a globally conserving solution as shown by the 0 residual line in Fig. 8.The MFA applied there is BC but the same result is obtained by any of the other algorithms.It also results in some reduction of the peak values of the field which is evident in Fig. 9.This can be explained if we consider that a MFA diagnoses that the total mass has been largely overestimated by cubic interpolation and has to remove mass to enforce conservation.As the mass is concentrated in a small area, across a few grid points, peak values will be inevitably reduced when the MFA is applied.Large interpolation errors as a result of large gradients and insufficient resolution near the source is the main reason for this mass overestimation.The sensitivity with respect to the specific mass fixer or quasi-monotone filter used was relatively small and all algorithms tested behave in a similar way.The biggest difference was found between the MG fixer and the remaining ones and this shows in Fig. 9.
Although it is difficult to obtain accurate results in test cases of advection of small-scale point sources with coarse (global)-resolution semi-Lagrangian models, useful qualitative results can still be obtained.The MFA may reduce the amplitude of the field but it will correct its total mass which is necessary for emission parameter estimation.

Conclusions
A MFA is a technique to correct the global mass conservation error that a non-formally conserving advection scheme introduces.It acts a posteriori to correct the solution after the field has been advected.In the context of a semi-Lagrangian scheme this means to correct the field after it has been interpolated to the departure point and before other source terms due to physical processes are added.
Different MFAs have been implemented (cf.Sect.3) in IFS based on different strategies for correcting the global mass conservation error.They all follow a weighted approach, i.e. weights are computed which determine how much to adjust each grid-point value.The aim is to correct the advected field in regions where the interpolation error is large.Results show that indeed these methods act in areas of steep gradients where the solution is not smooth while they apply very small corrections elsewhere.They achieve globally mass-conserving solutions without deteriorating accuracy at large scales.This has been demonstrated here by a set of 12-month forecast tests verified against ERA-Interim and standard 10-day forecasts at T511 L137 resolution verified against ECMWF operational analysis.A small local degradation of existing biases cannot be completely ruled out since the sign of the global mass error determines the sign of the corrections everywhere.The key results from this work are the following: 1.No significant differences have been found between the approaches at the hydrostatic scales tested.But there are small differences in cost.
2. Global conservation is achieved without deteriorating the solution.An exception is the volcanic plume case in which peak values are reduced.However, this sideeffect is also related to the lack of sufficient resolution.Despite this, global mass conservation is important for emission parametric estimates because the mass conservation error can reach up to half of the emitted mass.
3. The impact on forecast skill was neutral.
4. Noticeable impact was found from the type of quasimonotone limiter applied.In long integrations BS improves on the standard quasi-monotone scheme used in IFS.
Based on the above findings the recommendations on the use of the newly implemented MFAs in IFS are as follows: 1.For quasi-monotone cubic advection of moist quantities, BC is the preferred option as it is shape preserving and one of the cheapest.
2. If quasi-monotonicity is not essential and positivedefiniteness is sufficient, the cheapest fixer MG is sufficient.It is also the only one that can be applied for advection with linear interpolation and would be recommended for any model using such mixed approach.
3. The ZE fixer results in an accurate advection scheme and generates small increments.If quasi-monotonicity is not essential, it should be the best option for fields having background values away from zero.
4. Currently the BC fixer is recommended for simulations with chemical tracers because it is one of the cheapest and performs well in advecting correlated tracers (cf.Fig. 7).
5. For volcanic plumes, BC is also sufficient.
MFAs may be inappropriate at non-hydrostatic, cloudresolving scales.Future tests will include these regimes.
Ongoing developments in the PantaRhei project (ECMWF, 2013) will provide opportunities towards a strictly massconserving scheme for these regimes.Until such developments materialise, MFAs can provide a practical alternative for the applications supported by IFS and are attractive due to their low computational cost.

Fig. 1 .Figure 1 .
Fig. 1.Mass conservation error of the IFS SL advection scheme as a percent of initial global mass for ozone, Q, CLWC, CIWC at different horizontal and vertical resolutions using a quasi-monotonic bi-cubic or a linear (LIN, CLWC and CLIC only) interpolation scheme.

Fig. 2 .Figure 2 .
Fig. 2. Mass conservation errors as a percentage of initial global mass for Q, CLWC at T1279 L137 resolu forecast with/without PR and BC MFAs.

Figure 3 .
Figure 3. (a) Q and (b) BC fixer increment for Q (in kg kg −1 ) at t+24 h and 700 hPa height from a T1279 L137 forecast.

Figure 5 .
Fig. 5. 48 hrs timeseries of global rms-norms (left) and max-norms (right) of MFAs increments for CLWC expressed as a percentage of the rms-norm of the field.
appears) would indicate an improvement while cooling would indicate further deterioration.

Fig. 6 .Figure 6 .
Fig. 6.Experiments with BS limiter described in Table 1.Difference of vertical cross-sections of zonally averaged annual mean temperature fields.Plot (a): difference (in Kelvin) of control forecast from ERA-Interim.Plots (b-f): difference (in Kelvin) of experiments from control forecast.

Figure 8 .
Figure 8. Relative mass residual in volcanic plume simulations (SO 2 ) for different schemes Fig.9.Comparison of volcanic plume simulation with and without mass fixer using quasi-monotone cubic Lagrange at T1279 L91 resolution.The plotted quantity is the total SO2 content (in kg/m 2 ) per model gridpoint column.Experiments defined as in Table1.

Table 1 .
List of 12-month forecast experiments.Cubic qm setup using BS limiter instead of DEF.[unfiltered cubic, PR] Pure cubic Lagrange for moist fields, quasi-monotone advection by PR algorithm on moist fields.[cubic BSqm, BC] Cubic BSqm setup adding BC fixer on moist fields.[cubic BSqm, MG] Cubic BSqm setup adding MG fixer on moist fields.[cubic BSqm, ZE] Cubic BSqm setup adding ZE fixer on moist fields.