A New Mass Flux Correction Procedure for Vertically Integrated Energy Transport by Constraining Mass, Energy, and Water Budgets

Reconstruction of the atmospheric circulation from the numerical output of general circulation model (GCM) and reanalysis products often suffers from spurious values resulting in imbalances in the mass, energy, and water transport. A correction hence is needed prior to the analysis of energy transport. Most of the previous correction methods only employ a mass budget adjustment since it is the primary source of error. The present study proposes a new method of mass flux correction by constraining the mass, energy, and water budgets that can be applied to meridional energy transport and overturning circulation. The new method seeks a vertical dependent solution, unlike the conventional barotropic correction procedures. The procedure has been tested for one aquaplanet GCM and a reanalysis product. The method successfully corrects the spurious imbalances in the mean meridional energy transport and circulation.


Introduction
Meridional energy transport by the atmosphere-ocean system moderates the climate by cooling the hot tropics and warming the cold polar region. Hence, the spatial pattern of climate change and variability largely depends on the characteristics of heat transport. Although being such an essential diagnostic, the estimation of energy transport from general circulation model (GCM) outputs and reanalysis products suffers from spurious energy transport (e.g., Hill et al., 2017;Peters et al., 2008;Trenberth, 1991. A flux correction is generally required prior to the estimation. Since the budget imbalances in the mass, energy, and water budgets are important sources of errors in the energy transport estimation, it is highly desirable to formulate a standard flux correction procedure that satisfies these budgets. The poleward energy transport is usually estimated from either direct or residual perspective. In a direct estimate, the poleward energy transport is calculated as the amount of energy carried by air parcels. In contrast, in the energy budget-based residual estimate, the transport of energy is calculated by inferring the required divergence of energy at each latitudinal belt to meet the energy balance of the top of the atmosphere (TOA) and surface. Both the direct and residual methods encounter imbalances in the estimation of energy transport. The direct method encounters unexpected estimates of the value and structure of heat transport that is inconsistent with the column energy balance. Contrastingly, the residual estimation often ends up with The spurious imbalances in the mass, energy, and tracer budgets in the outputs of GCMs and reanalysis stem primarily from the interpolation of model output from their native coordinates to regular latitude-longitude and pressure coordinates, making the mass continuity no longer satisfied (e.g., Trenberth, 1991. The model's poor vertical resolution is another source of error, especially in the bottom boundary layer and in the upper troposphere. Although the model equations conserve mass and energy budgets, the numerical model outputs often do not conserve them even in the native grids. This imbalance usually arises from the many approximations and discretization techniques used in various numerical codes representing diffusion, advection, spatial time filters, and time differencing (e.g., Becker, 2003;Jablonowski & Williamson, 2011;Lucarini & Ragone, 2011). In reanalyses, imbalances are also induced by the assimilation of uneven observation and low-frequency sampling of outputs. The imbalances in individual conserved quantities do not merely affect the respective budgets, but they distort other dependent budgets and many essential climate diagnostics that involve mass weighting and are based on residuals of budget equations (e.g., Boer, 1982;Fortelius & Holopainen, 1990;Trenberth, 1991).
There have been many ad hoc correction procedures to fix the issue of imbalances in the budgets and the consequent spurious energy transport (e.g., Boer, 1982;Fortelius & Holopainen, 1990;Lembo et al., 2019;Peters et al., 2008;Trenberth, 1991;Trenberth et al., 2001). Since a greater portion of the imbalances in energy transport is associated with the mass budget imbalance, the first-order correction to energy budget is to adjust the mass budget. Such a correction removes most of the imbalances in the energy transport. However, correcting the mass budget alone does not guarantee a reasonable correction in energy transport in many instances since the transport also depends on the energy and water budget.
Here, we propose a new method of flux correction that is applicable for both the zonal mean meridional energy transport and circulation by constraining the mass, energy, and water budget. Unlike most of the previous correction procedures, where a barotropic correction is applied at all levels, the new method seeks a vertically dependent correction by using three different vertical scalings: pressure scale height, equivalent depth of tropospheric baroclinic mode, and water vapor scale height. The method combines both the direct and residual methods of energy transport estimation to formulate the correction procedure. The direct estimation method requires very high-frequency instantaneous energy and velocity fields to compute the energy fluxes. In contrast, the residual estimation method does not require high-resolution instantaneous data since the surface and TOA energy balance terms involved in the residual method are accumulated quantities. However, the residual method cannot decompose the transport into mean and eddy parts, unlike the direct estimation. By combining both direct and inferred transport and by constraining the mass, energy, and water conservation, the new procedure enables direct estimation of flux corrected meridional heat transport and its various components. The mass flux correction can also be applied to the estimation of the mean meridional overturning circulation.

Correction Procedure
The proposed methodology seeks to find a correction to the meridional velocity field to ensure mass, energy, and water vapor conservation. We assume the correction velocity (v c ), an exponential polynomial form, as where v c is a function of both pressure (p) and latitude ( ). v 1 , v 2 , and v 3 are the coefficients which have only latitudinal dependence. Z * is a height function defined as where H 0 = RT 0 ∕g ≈ 8-10 km is the pressure scale height of an isothermal atmosphere. We choose three different heights corresponding to the pressure scale height (H 1 ≈ 10 km), the equivalent depth of tropospheric baroclinic mode ( H 2 ≈ 5 km) and water vapor scale height (H 3 ≈ 2.5 km). These scales, H 1 , H 2 , and H 3 , are chosen to approximately scale the correction corresponding to the scales associated with the transport of mass, dry static energy, and water vapor, respectively. All these three exponentially decaying scales aloft helps to limit the correction in the stratosphere, where the mean circulation should be weak compared to the troposphere. Combining Equations 1 and 2, the expression for v c can be rewritten as The exponents now are 1 = We can now find v 1 ( ), v 2 ( ), and v 3 ( ) such that the corrected velocity field has no net mass flux and match the energy and water vapor transport from the energy and water budget.
Column integrated zonally averaged budgets for mass, energy, and water vapor at equilibrium are respectively. The overbars and the square brackets denote the time mean and the zonal mean, respectively. Here, a is the Earth's radius, g is gravitational acceleration, v is the meridional velocity, m is moist static energy (MSE), and q is the specific humidity. The MSE is given by m = C pd T + gz + L v q, where the first, second, and third terms are the dry-air enthalpy, gravitational potential energy, and latent heating, respectively. The dry air enthalpy and the gravitational potential energy together constitute the dry static energy (DSE). Here, C pd is the dry air heat capacity at constant pressure, T is temperature, z is height, and L v is the latent heat of vaporization. In Equation 5, we neglect the transport of kinetic energy. The terms F m and F q are the zonally averaged column energy and water transport inferred from (residual estimation) the column energy and water budget, respectively. The annual mean budget residuals are provided in the supporting information (Text S5).
The residual estimate for the energy transport F m can be derived from the surface and TOA energy budget as follows, where E is the zonally averaged column integrated heat content and  net is the net forcing within the atmosphere, that is, (R toa + R s c +  lh +  sh ). Here, R toa and R sfc are the net radiation at TOA and at the surface, respectively. C vd is the specific heat of dry air at constant volume. And  lh and  sh are the surface latent and sensible heat fluxes, respectively.
Since the energy budget retrieved from the model has an imperfect closure, the net heat transport derived from the radiation balance may not converge to 0 at the poles. Therefore, we need to correct the energy budget imbalance prior to the estimation of inferred energy transport. To correct the net energy out to space, "Carissimo correction" (Carissimo et al., 1985;Lucarini & Ragone, 2011;Lucarini et al., 2014) is applied by assuming the imbalance is uniform throughout the latitude. We first average the zonally averaged net radiation in the atmosphere throughout the latitude, which in principle is the radiation imbalance. The imbalance is then subtracted from the zonally averaged radiation budget at each latitude. It is worth noting that the Carissimo correction is sufficient for the estimation of inferred MSE transport by total circulation, but it does not allow the decomposition of transport by different circulation components. Hence, we need a correction procedure for the direct estimation of transport, where the transport can be decomposed into the mean and eddy part.
Similar to the inferred MSE transport, the zonally averaged column water transport can also be inferred from the column water balance as follows, Here, the source and sink terms for the water budget,  and , are the rate of evaporation and precipitation, respectively.q is the vertically integrated specific humidity.
At equilibrium, the energy and water vapor tendency terms in Equations 7 and 8, respectively, become 0. This is mostly true in an annual average scenario, but seasonally the energy tendency term may not be negligible due to the heat capacity of the ocean. So, the tendency term should be involved while estimating the seasonal heat transport.
Substitute F m and F q back to Equations 5 and 6, respectively. Now, Equations 4-6 constitute a linear algebraic system with three equations and three unknowns (v 1 , v 2 , and v 3 ), which can be solved at each latitude to get the velocity correction, v c (see supporting information, Text S1.).

Data
We choose an aquaplanet GCM and a reanalysis to test the validity of the procedure. Aquaplanet simulation has been chosen to exploit its simple structure to better understand the influence of the procedure. Unlike full-fledged GCM and reanalysis, the aquaplanet simulations are hemispherically symmetric and do not have complexity related to unequal land distribution. Such a framework helps to understand the effect of correction easily. However, to prove the functionality of the procedure in full-fledged operational data sets, we have tested the procedure in the ERA-Interim Reanalysis. Three years (2000)(2001)(2002)(2003) of data were analyzed. Six-hourly data have been used for the direct estimation of energy transport in both GCM and reanalysis.
An aquaplanet version of High Resolution Atmospheric Model (HiRAM) (Anderson & Coauthors, 2004;Zhao et al., 2009), developed at Geophysical Fluid Dynamics Laboratory, has been used to conduct experiments. HiRAM is built by modifying the atmospheric component of GFDL GCM, Atmospheric Model version 2 (AM2). Major updates in HiRAM are the use of a quasi-uniform horizontal grid spacing in finite-volume core using a cubed-sphere grid topology (Putman & Lin, 2007) instead of the latitude-longitude grid of finite-volume dynamical core (Lin, 2004) in AM2, improved resolution, and simplified parameterizations for moist convection and large-scale (stratiform) cloudiness. The present simulations are run at a resolution of ≈50 km, and the output is interpolated into regular latitude-longitude and pressure coordinate. HiRAM is unique in its moist convection. Following the idea of resolving the convection with the help of parameterized convection in high-resolution simulations, proposed by Pauluis et al. (2006) and Garner et al. (2007), HiRAM uses a scheme based on a shallow convection parameterization by Bretherton et al. (2004) instead of the relaxed Arakawa-Schubert convective closure (Moorthi & Suarez, 1992) in AM2. Such an approach makes the parameterized convection less intrusive and allows explicit convection to take over the vertical flux. A more detailed description of the model can be found in Zhao et al. (2009).
The current aquaplanet version of HiRAM uses a simple mixed-layer ocean as the bottom boundary (Kang et al., 2008). The mixed layer is a static slab ocean with a specified heat capacity, which defines the depth of the mixed layer. The present study uses a heat capacity of 8 × 10 7 J m −2 K −1 , corresponding to a depth of ≈20 m. Small heat capacity or depth has been chosen to attain equilibrium quickly. The model employs a seasonal cycle with incoming daily mean solar radiation computed as a function of latitude and day of the year. There is no diurnal cycle in the insolation. We have conducted aquaplanet simulations for 13 years. Only the last 3 years were taken for the analysis to make sure the analysis is done in an equilibrium climate.

Results
In order to examine the effectiveness of the correction procedure, zonally averaged column integrated MSE, DSE, and latent heat (LH) transport by total circulation (both direct and residual method ) estimated before and after the mass flux correction are portrayed (Figure 1). The left column (Figures 1a-1c) represents the estimates from HiRAM aquaplanet simulations, while the right column (Figures 1d-1f) represents the estimate from ERA-Interim reanalysis. In HiRAM simulations, MSE transport computed directly from the model output (Figure 1a, solid green line) has a double peak structure in both hemispheres: one at ∼20 • and the other at ∼45 • latitude overestimating the transport. However, such a double peak structure is inconsistent with the energy balance at TOA. The energy balance dictates a seamless transport with a single peak in the midlatitude (∼35 • latitude), as is evidenced by the transport inferred from the column energy balance (blue dashed line). The mass flux correction procedure efficiently rectifies the disparity in the direct estimation of transport (solid red line). The correction helps to get rid of the physically inconsistent double-peak structure.

Geophysical Research Letters
10.1029/2020GL089764 Figure 1d shows MSE transport in ERA-Interim reanalyses, before and after correction. Unlike the HiRAM simulation, the bias in the direct estimation of MSE in reanalysis has an interhemispheric asymmetry. The direct estimation before correction (solid green line) underestimates the transport in the Southern Hemisphere except near the South Pole. In contrast, the direct estimation overestimates the total energy transport in the Northern Hemisphere subtropics and midlatitude. Therefore, the correction involves transport of energy from the Northern Hemisphere to the Southern Hemisphere. Moreover, the difference between directly computed transport before the correction and the corrected inferred transport is higher in reanalysis compared to HiRAM simulations. A key reason for this disparity is the fact that the aquaplanet simulation is devoid of important error sources such as complex topography and data assimilation. However, irrespective of the complexities, the correction procedure effectively corrects the direct estimated MSE (solid red line) to agree with the required transport (dashed blue line). It has also been shown ( Figure S1 and Text S2; see supporting information) that the new procedure has significantly improved the energy transport estimation compared to the conventional barotropic method. The new correction procedure is also effective in the seasonal estimate of the transport (Figures S1 and S2 and Text S2; see supporting information).
By construction, the Carissimo correction fixes the unphysical residual energy transport out to space in the residual transport (dashed magenta and blue line) first. Later the velocity correction leads to an exact match between the (corrected) direct and residual transport for both aquaplanet simulations and reanalysis. It is worth understanding how the procedure corrects DSE and LH components of the total energy transport. Overall, the correction mainly happens in the DSE component of heat transport, while the correction in LH total transport is small. Therefore, the pattern of correction in total transport of MSE is dominated by the correction in the DSE transport.
In the HiRAM case, the DSE total transport (Figure 1b) computed directly from the model output overestimates the transport in the tropics and extratropics. The influence of correction on the total LH transport ( Figure 1c) is mainly over the midlatitude, but the magnitude of correction is still very small compared to DSE correction. The correction procedure effectively fixes the imbalances in both DSE and LH total transport.
The direct estimate of the ERA-Interim DSE total transport (Figure 1e) underestimate the transport in most of the Southern Hemisphere except near the Antarctic. Near the Antarctic, direct estimation overestimates DSE transport. In the Northern Hemisphere, the bias is relatively less compared to the Southern Hemisphere and is confined to the tropics. The correction in the LH transport ( Figure 1d) is relatively less. The direct estimate of LH transport overestimates the poleward transport at both midlatitudes. However, the error in the equatorward transport at the tropics is quite small.
Next, we examine how the correction procedure works in the mean and eddy components of the MSE, DSE, and LH transports. To perform this analysis, first, the velocity and energy fields are decomposed into time mean and transient eddy components. The transient eddy components are derived by subtracting the time mean component from the total field. The correction mainly occurs in the Eulerian mean component of the transport, while it is negligible in the eddy part ( Figure S1 in the supporting information). Figure 2 displays the MSE, DSE, and LH transport by the Eulerian mean circulation. Figures in the left column (Figures 2a-2c) are the Eulerian mean energy transport components from HiRAM simulations. In contrast, those in the right column (Figures 2d-2f) are the ERA-Interim counterpart.
In the case of Eulerian mean MSE transport in HiRAM aquaplanet simulation (Figure 2a), the correction reduces the directly estimated poleward MSE transport associated with the Hadley cell. However, the procedure enhances the equatorward transport by Ferrel cell mostly by broadening the cell poleward, which in turn contributes to weaken the extratropical poleward MSE transport by the total circulation (Figure 1a). The decomposition of the MSE transport into DSE and LH components shows that the correction in mean MSE transport is dominated by the correction in the DSE component (Figure 2b). The pattern is very similar to that of MSE transport. That is, the correction reduces the poleward transport associated with the Hadley cell and increases the Ferrel cell energy transport. The correction in the Eulerian mean LH transport (Figure 2c) is small compared to that of DSE. Since the correction in the energy transport by eddies and the LH transport by Eulerian mean circulation are negligible, the correction in the MSE total transport is dominated by the correction in the Eulerian mean DSE transport (Figure 2b).

Geophysical Research Letters
10.1029/2020GL089764 Figure 2. Same as Figure 1, but for transport of energy by Eulerian mean circulation.
The correction in the Eulerian mean energy transport in the ERA-Interim case (Figures 2d-2f) is different from the HiRAM case. The correction is maximum in the southern hemispheric tropics. Unlike that of the aquaplanet case, the correction is to enhance the poleward transport of DSE and MSE by the southern hemispheric Hadley circulation. Contrastingly, correction weakens the MSE and DSE transport in the northern hemispheric tropics. In the midlatitudes, the correction weakens the equatorward transport of DSE and MSE by the southern hemispheric Ferrel cell. The correction is comparatively smaller in the northern hemispheric midlatitude. The correction weakens the equatorward LH transport (Figure 2e) by the southern hemispheric Hadley cell and enhances the poleward LH transport by the southern hemispheric Ferrel cell. In the northern hemispheric midlatitude, the correction enhances the poleward LH transport. There is To understand how the velocity correction fixes the energy transport anomalies, correction velocity (v c ) obtained by solving the linear system of equations (Equations 4-6) and the original meridional velocity (v) are shown in Figure 3. The correction coefficients v 1 , v 2 , and v 3 are provided in the supporting information ( Figure S2). An experiment of velocity correction sensitivity to the height scales' choices demonstrates ( Figure S6 and Text S4; see supporting information) that the proposed physically motivated heights (H 1 ≈ 10 km, H 2 ≈ 5 km, and H 3 ≈ 2.5 km) are the optimal choice. In the HiRAM case (Figures 3a and 3b), the correction acts to make the tropical low level equatorward flow more shallow, while weakening the overturning Hadley circulation. As such, the correction decreases the poleward transport of DSE while leaving the moisture transport mostly unaffected. In the midlatitudes, both the low-level poleward flow and the equatorward upper-level flow are enhanced by the correction, which strengthens the Ferrel cell and expands it poleward. Such a circulation change is consistent with the enhanced equatorward DSE transport by Eulerian mean circulation ( Figure 2b) and a weakening of the poleward MSE transport by total circulation (Figure 1a).
Unlike the HIRAM case, the correction in ERA-Interim (Figures 3c and 3d) is hemispherically asymmetric. The correction acts to weaken the northern hemispheric Hadley cell and enhances the southern hemispheric counterpart. In other words, the correction will decrease the mean poleward DSE transport in the Northern Hemisphere and vice versa in the Southern Hemisphere. In the midlatitude, the correction weakens the Ferrel cell and extends poleward in the Southern Hemisphere consistent with the weakened Eulerian mean DSE transport (Figure 2e) and the enhanced poleward transport by total circulation (Figure 1d).
It has been demonstrated that the proposed mass flux correction method works efficiently with the meridional energy transport and its dynamic and energetic components. However, it is also essential to assess how the velocity correction (Figure 3) impacts the mean meridional circulation. In other words, we have to validate the potential of the present correction procedure for the mean meridional circulation. The impact of the velocity correction on the mean meridional circulation in both isobaric and isentropic coordinates are portrayed in Figure 4. Only the results from the HiRAM aquaplanet simulation is provided for brevity. In isobaric coordinates, the mass correction reduces the strength of tropical circulation and increases that of the extratropical cell (Figures 4a and 4b). The correction also broadens the Ferrel cell poleward. Vertically, the correction mostly confined to the troposphere. In the tropics, the correction peaks in the middle to upper troposphere, while the correction is mostly at the lower troposphere in the midlatitude. The correction is about 10% of the stream function in the peak locations.
In isentropic coordinate (potential temperature coordinate) (Figures 4c and 4d) also, the velocity correction reduces the strength of tropical circulation and increases the strength of extratropical circulation. The correction also broadens the extratropical cell further poleward. Also, the correction outside the tropics peaks in the subtropical upper troposphere in the isentropic coordinate. The correction weakens the large-scale subsidence over the middle to upper subtropical troposphere. The midlatitude circulation is enhanced in both hemispheres as a consequence of the correction. Specifically, the surface branch of poleward mass transport from the midlatitude is enhanced. Although the correction in the tropics is about 10% of the total mass transport as in the isobaric coordinate, in the lower stratosphere, the correction is comparable to the original stream function. This is likely due to having the corrected velocity to be fairly deep (with H 1 ≈ 8 km) thus allowing it to extend in the to stratosphere.
In general, irrespective of the coordinate system (isobaric or isentropic coordinates), the correction in the mean meridional circulation is consistent with the energy transport correction, which was not possible in the conventional barotropic correction. The weakening of the tropical overturning and the strengthening of the extratropical circulation is consistent with the Eulerian mean DSE changes (Figure 2b), which is the dominant correction in the total meridional energy transport. However, while the correction ensures that the tropospheric circulation is consistent with the energy and water budget, it may negatively impact the estimate of the stratospheric overturning.

Conclusions
The present study proposed and validated a new post-hoc mass flux correction procedure specifically to apply for the column integrated zonally averaged meridional energy transport and the mean meridional overturning circulation. Unlike the conventional methods, where the correction is employed by constraining the mass conservation, the present method takes into account the energy and water conservation in addition to the mass conservation. Moreover, the present method implements a vertically dependent solution instead of the conventional barotropic solution.
The proposed method of correction has been validated using the simulation from an aquaplanet general circulation model and reanalysis product. It has been shown that the method effectively corrects the spurious imbalances in meridional energy transport. Moreover, the method satisfies the major budget constraint in the climate. It has also been demonstrated that the method is equally useful in correcting the mean meridional overturning circulation in both isobaric and isentropic coordinates. The procedure can potentially be applied to reanalyses and GCM outputs.