Parameterizing the vertical downward dispersion of ship exhaust gas in the near field

Estimating the impact of ship emissions on local air quality is a topic of high relevance, especially in large harbor cities. For chemistry-transport modeling studies, the initial plume rise and dispersion play a crucial role for the distribution of pollutants into vertical model layers. This study aims at parameterizing the vertical downward dispersion in the near field of a prototype cruise ship, depending on several meteorological and technical input parameters. By using the microscale chemistry, transport and stream model (MITRAS), a parameterization scheme was developed to calculate the downward dispersion, i.e., the fraction of emissions, which will be dispersed below stack height. This represents the local concentration in the vicinity of the ship. Cases with and without considering the obstacle effect of the ship have been compared. Wind speed and ship size were found to be the strongest factors influencing the downward dispersion, which can reach values up to 55 % at high wind speed and lateral wind. This compares to 31 % in the case where the obstacle effect was not considered and shows the importance of obstacle effects when assessing the ground-level pollution situation in ports.

1 Introduction

Ship emissions are among the harmful anthropogenic influences on air quality and human health, especially in big harbor cities. Regarding air quality, this covers various gaseous pollutants like sulfur dioxide (SO₂) and oxides of nitrogen (NO_x), as well as particulate matter (PM). Corbett et al. (2007) presented a study on the global effect of ship-related particulate matter emissions on human health and found that they are at least partially responsible for around 60 000 cardiopulmonary and lung cancer deaths annually, mainly in coastal regions of Europe, and south and east Asia. Despite numerous measures to reduce these impacts, like the International Maritime Organization's global sulfur cap to 0.5 % maximum sulfur content in marine fuels since 1 January 2020 (IMO, 2021), research on the impacts of ship emissions on human health and air quality remains an ongoing topic (Barregard et al., 2019; Broome et al., 2016; Lin et al., 2018; Liu et al., 2016; Ramacher et al., 2019; Sofiev et al., 2018; Zhang et al., 2019).

The urge to quantify ship emissions gained even more attention due to the fast growth in shipping activity during the last decades (Brandt et al., 2013; United Nations Conference on Trade and Development, 2019). Many ship plume modeling studies focus on global- or regional-scale plume dispersion, chemistry and their parametrization (e.g., Aksoyoglu et al., 2016; Huszar et al., 2010; Vinken et al., 2011). For emission calculation, they can make use of global (Corbett et al., 2007; Wang et al., 2008) or regional (Aulinger et al., 2016; Jalkanen et al., 2009, 2012) shipping emission inventories. Increasing trends in ship emissions in northern Europe have been modeled for the North Sea (Matthias et al., 2016) and Baltic Sea (Karl et al., 2019), and the efficiency of emission reduction measures has been evaluated.

Matthias et al. (2018) pointed out the importance of correct spatial and temporal distribution of emissions in chemistry-transport models, which are connected with uncertainties. The distribution depends strongly on data availability, interpolation procedures and initial assumptions. For example, an emission overestimation for single ships can occur if ship emissions are diluted instantaneously and equally into a large grid (von Glasow et al., 2003; Vinken et al., 2011). This problem can be overcome by using smaller-scale models and a bottom-up approach for emission inventories (Aulinger et al., 2016; Eyring et al., 2010; Jalkanen et al., 2009, 2012; Johansson et al., 2017).

However, new problems arise when modeling is performed on a single-ship level. On the technical side, this includes the correct localization of ships, which can nowadays be done by using automatic identification systems (AISs). The emission factors are calculated as a function of several technical parameters like fuel type, engine load, engine power, engine type and ship size, ship class, as well as the performed activity (e.g., cruising, maneuvering, or berthing), the latter to be derived from AIS data (Aulinger et al., 2016).

But also new methods for spatially allocating the emissions from big ships need to be developed, since the emission height often lies significantly above ground level. The exhaust gas leaves the stack with a certain exit velocity and a temperature of several hundred ^∘C. These quantities depend on technical parameters of the individual ship, which are often unknown.

Many analytical single-plume models are based on Gaussian dispersion formulas (Briggs, 1982; Janicke and Janicke, 2001; Schatzmann, 1979). This means that the pollutant distribution corresponds to a normal probability distribution. An example is the offshore and coastal dispersion (OCD) algorithm of Hanna et al. (1985). An accurate representation of plume rise and downward dispersion processes in the near field under different meteorological conditions is important, since it changes the effective emission height and may cause the vertical concentration profile to deviate from a Gaussian shape (Bieser et al., 2011; Brunner et al., 2019). Based on a large eddy simulation study, Chosson et al. (2008) pointed out that Gaussian plume dispersion models might not be well suited for the early plume development.

Despite not running at full engine power inside of the harbor, ocean-going ships still consume large amounts of fuel for heat and electricity production and therefore emit atmospheric pollutants while at berth (Hulskotte and Denier van der Gon, 2010). These have been found to be up to 5 times higher compared to other activities like maneuvering or cruising during the course of a year, as ships spend more time at berth and also have a high auxiliary engine power demand for hotel services (Tzannatos, 2010). This can lead to severe air quality problems in harbor areas. Murena et al. (2018) applied a computational fluid dynamics model to assess the impact of cruise ship emissions on the facades of waterfront buildings in Naples, Italy. The highest SO₂ concentrations were modeled for areas near the ships and on the port-front facades of the first line of buildings.

More research is needed to better understand the effect of technical and meteorological parameters on the downward dispersion process that causes these strong pollution scenarios inside of a harbor. Although the OCD model of Hanna et al. (1985) includes effects of pollutant downward dispersion behind the obstacle, i.e., the vessel, by lowering the effective plume height and adjusting dilution parameters in the model, this effect has yet to be applied to large ships.

The aim of this study is to quantify the main factors that have an impact on the downward dispersion process for a large cruise ship. For this purpose, an Eulerian microscale model is used to calculate the downward dispersion in close proximity to the ship. Furthermore, a parameterization is developed for the downward dispersion depending on the crucial meteorological influences and the technical specifications of the ship. Finally, it will be shown under which conditions the obstacle effect on the downward dispersion needs to be considered.

2 Methodology

The dispersion of an exhaust plume is affected by several meteorological and technical parameters (Fig. 1). The upward movement, i.e., the plume rise, is mainly determined by the initial temperature of the exhaust and its exit velocity, which can be calculated by dividing the gas volume flow by the stack diameter. The stack angle describes whether the exhaust flow is directed vertically, horizontally or at an angle. The stack height only has an indirect effect on the plume rise, as higher emitted gases experience a stronger wind speed inside the boundary layer.

Figure 1 Conceptual model of parameters affecting the shape and movement of a ship plume (bold text indicates technical parameters; italic indicates ambient parameters).

Turbulence enhances the plume dispersion, leading to dilution of the embedded gases by entrainment of ambient air into the plume. The dispersion increases with the wind speed. It depends also on the ship geometry and the flow direction of the wind towards the vessel. Furthermore, a stronger turbulence occurs in the case of higher surface roughness.

The ambient vertical temperature profile determines the atmospheric stability. The presence of an inversion can strongly decrease the strength of the plume dispersion, as it thermodynamically hinders the vertical movement of air masses. Depending on the altitude of the inversion and the exhaust temperature, the plume may or may not break through the inversion.

2.1 MITRAS

The microscale chemistry, transport and stream model (MITRAS) is a non-hydrostatic, three-dimensional Eulerian model, based on the Navier–Stokes equations, the continuity equation and the conservation equations for scalar properties like temperature, humidity and trace gas concentrations (Grawe et al., 2013; Salim et al., 2018; Schlünzen et al., 2003, 2018). It accounts for obstacle-induced turbulence on the wind field as well as effects of thermal stratification.

In this study, a non-equidistant grid is used with the highest resolution of 2 m × 2 m × 2 m close to the ship. The chosen domain has an overall size of around 1 km × 1 km horizontally and 500 m vertically. The surface cover for the whole domain is water and the roughness length is calculated from the wind speed (see Schlünzen et al., 2018, for detailed equations). Its values are close to zero.

The emission occurs continuously in one model cell right above the ship stack, which is an impenetrable obstacle cell (Fig. 2). The emitted gas is as a passive trace gas (e.g., CO₂ or non-reactive SO₂). No chemical reactions occur in the simulations. The emission cell has a constant temperature, which corresponds to a given exhaust temperature and a vertically directed exhaust velocity. The wind field is affected by Coriolis force and friction force, which cause the wind to slightly turn counterclockwise according to an Ekman spiral. Furthermore, the flow field is modified by the obstacle itself, the high temperature of the exhaust and the exit velocity. No deposition occurs in the model domain, the surface is a mirror source which reflects the concentration when the lowest model layer is reached.

Figure 2 Visualization of the stack emission for wind direction from left to right. Passive trace gas emission occurs in the cell above the stack, which has a constant exhaust temperature and a vertically directed exhaust velocity. The arrows indicate the change of the ambient wind field due to the obstacle and the plume temperature.

2.2 Meteorological data

Idealized meteorological conditions are used to investigate effects of single variations of input parameters on the dispersion process. The range of input values is listed in Table 1. One input parameter per model run was varied while the other meteorological and technical parameters were fixed at predefined default values.

Table 1 Input parameters for this study. While varying a single input parameter in the investigation range, all others remain at default setting.

The ambient temperature is set to 15 ^∘C at the surface. It changes with altitude according to the given ambient temperature gradient, which represents the atmospheric stability. The value of ambient temperature itself has a negligible effect on the plume dispersion compared to the plume temperature and was therefore not varied in this study.

The atmospheric stability is varied in a range of different lapse rates, covering one unstable condition (−1.2 K × 100 m⁻¹), one neutral condition (−0.98 K × 100 m⁻¹) and several stable conditions including inversions (up to +0.5 K × 100 m⁻¹).

The wind speed is investigated in the range of 2–15 m s⁻¹. The limits were chosen according to hourly wind speed data from the Hamburg weather mast in 2018 (see Appendix A) and can also be seen as representative for other large northern European ports including Rotterdam and Antwerp. The value 2 m s⁻¹ is close to the 5th percentile and 15 m s⁻¹ corresponds to the 95th percentile at 280 m measurement height. This covers most of the naturally occurring scenarios. The selected default value is 5 m s⁻¹ which fits well with the mean wind speed in Hamburg at a height of 50 m, which is close to the stack height.

The effect of wind direction is relevant in correspondence to the orientation of the ship. Frontal wind is herein defined at an angle of 0^∘ and lateral wind at 90^∘. Oblique wind conditions lie between these values.

2.3 Ship characteristics

This study represents a cruise ship prototype. From an online database (Port of Hamburg, 2020), the average length and width of cruise ships that were visiting Hamburg harbor during the years 2018 to 2019 has been calculated. The stack height was approximated from freely available photos (e.g., Vesseltracker, 2020). The ship prototype has a length of 246 m, a width of 30 m and a stack height of 52 m (see Table 1 and Fig. 3). This corresponds to a typical cruise ship that can carry between 1000 and 2500 passengers. A non-moving source is assumed, i.e., a hoteling ship at berth.

Figure 3 Side view of the prototype cruise ship in the MITRAS domain with the x axis located at the stack position.

The study goes beyond a case study. A loaded container ship of similar size, and exhaust characteristics would deliver similar results because its shape is comparable. On top of that, for all investigated input characteristics, the results of stack-only cases are presented as well. Therefore, one can assume that results for smaller ships lie between these two cases.

The exhaust gas temperature depends on technical parameters of the ship's engine and can be found in engine data sheets provided by manufacturers like Caterpillar (CAT, 2020), Wärtsilä (Wärtsilä, 2020) and MAN (MAN, 2020) on their websites. For large cruise vessels, it ranges between approximately 300 and 400 ^∘C, depending on the used engine power. However, the exhaust temperature can be lowered by 75–100 ^∘C when a heat exchanger which generates electric energy from the excess heat is in operation (Murphy et al., 2009). Therefore, the temperature effects are investigated for 200, 300 and 400 ^∘C plumes to cover a realistic spectrum. Similarly, the exit velocity was assumed from these data sheets. It depends on the engine type (main engine or auxiliary engine) and the used engine power and was investigated in a range of 4–12 m s⁻¹.

2.4 Plume dispersion in different regimes

When investigating plume dispersion, one needs to separate two regimes: the momentum-driven regime and the buoyancy-driven regime. In the momentum-driven regime, the movement of the plume is affected by (a) the initial plume rise due to both the exit velocity and the high-temperature convective upward transport and (b) the dispersion due to turbulence generated by the obstacle (i.e., the ship) inside the wind field. In the buoyancy-driven regime, the movement of the plume is determined by the wind field and turbulence generated by the ambient conditions (e.g., orography effects and surface roughness). Here, the plume temperature is equal to the ambient temperature. The microscale model (MITRAS) can investigate plume behavior in both regimes on a high resolution.

MITRAS is used to capture the initial plume rise and turbulence effects in the momentum-driven regime. The vertical concentration profiles are calculated at a distance outside of the momentum-driven regime, i.e., when the buoyancy-driven regime is reached. Then, the concentration profiles are calculated on a 100 m × 100 m area column with layer-mean values (Fig. 4). The calculation of these column values has two benefits. First, it covers the mean behavior of the whole plume better than single values of 2 m × 2 m × 2 m grid sizes, since the movement of the plume can be highly variable. Second, the concentration profiles can then also be transferred into larger-scale models which usually have a much coarser grid (e.g., 100 m × 100 m horizontally). However, the coupling of MITRAS results into a larger-scale model will be part of a future study and is not covered here.

Figure 4 Schematic sketches of the investigation area. Vertical concentration profiles are evaluated at a distance of 100 m downwind of the ship for layers of 100 m × 100 m.

Since the plume needs to have cooled down to ambient temperatures to be considered outside the momentum-driven regime, test simulations have been performed to find a distance at which this condition is met (see Appendix B). This was the case at a distance of 100 m downwind of the ship. Therefore, all concentration profiles are calculated as 100 m × 100 m columns with average concentration per layer at a distance of 100 m downwind of the ship.

In the following, the term “downward dispersion” (D) is defined as the relative proportion of the total concentration column in the layers below the stack height.

$$\begin{array}{}\text{(1)}& D={\displaystyle \frac{\underset{\mathrm{0}}{\overset{h_{\mathrm{stack}}}{\int }}c}{\underset{\mathrm{0}}{\overset{h_{\mathrm{top}}}{\int }}c}}\cdot \mathrm{100}\mathit{\%},\end{array}$$

where h_top is the altitude of the highest model layer (500 m), h_stack is the stack height (52 m), and c is the total concentration. A mean downward dispersion is calculated for the described 100 m × 100 m column at a distance of 100 m downwind from the stack. From an application perspective, this downward dispersion parameter is an indicator of the pollution situation in the vicinity of the ship and useful to evaluate the level of pollution inside of a harbor.

For single regression analyses, downward dispersion values are investigated depending on the variation of a single input parameter at a time while the others remain at default settings (Table 1). To assess the sensitivity of the downward dispersion to each input parameter, an effective range r is calculated. It is defined as the difference between the highest and the lowest downward dispersion values for one regression:

$$\begin{array}{}\text{(2)}& r_i=D(max{)}_i-D(min{)}_i,\end{array}$$

where i is the individual input parameter that is varied while the other remains at the default setting. The effective range describes how strongly one parameter can change the downward dispersion and helps to evaluate which input parameter has the strongest impact in the given range of values.

3 Results and discussion

Figure 5 presents an exemplary output of the MITRAS model for the default conditions, i.e., frontal wind at 5 m s⁻¹, exit velocity of 10 m s⁻¹, exhaust temperature of 300 ^∘C and an ambient temperature gradient of −0.65 K × 100 m⁻¹. The concentration values result from an emission of 50 kg trace gas per hour.

Figure 5 MITRAS model results for default conditions (frontal wind at 5 m s⁻¹, exit velocity 10 m s⁻¹, exhaust temperature of 300 ^∘C and an ambient temperature gradient of −0.65 K × 100 m⁻¹).

The following subsections describe the results of single- and multi-parameter regressions that were performed in order to describe the relationship between the downward dispersion and the input parameters. From the multi-parameter regression, a parameterization is derived that covers all input parameters in the investigation range. A bootstrapping procedure is presented to test how well the parameterization results match with the MITRAS model results. The obstacle effect is evaluated, and finally some limitations of the modeling approach are discussed.

3.1 Results of single-parameter regressions

Single-parameter regressions are performed after basic statistic formulae (see Appendix C) to investigate the impact of individual input parameters, i.e., wind speed, exit velocity, wind direction, plume temperature and atmospheric stability on the downward dispersion.

3.1.1 Effect of wind speed and exit velocity

The dependence of the downward dispersion from wind speed was modeled in the range of 2–15 m s⁻¹ at the uppermost model layer, which is set as the input parameter. It is slightly lower at stack height following the logarithmic vertical wind profile. See Table C1 for the exact values at stack height.

Figure 6 presents results of a single linear regression for the dependence of downward dispersion on varying wind speeds with and without the obstacle effect. Other input parameters remained constant at default values (Table 1). A linear relationship with correlation coefficients R² of 0.98 was found for both runs with and without the ship, respectively. At high wind speeds, the turbulence behind the obstacle causes strong downward dispersion. Under these settings, the wind speed has an effective range on the downward dispersion of 40.3 % with and 21.1 % without the ship, making the wind speed a crucial factor influencing the downward dispersion (Fig. 6 and Table 2).

Figure 6 Dependence of the downward dispersion D on different wind speeds v_wind with and without the obstacle effect.

Table 2 Effective ranges of investigated input parameters on the downward dispersion under default settings.

A similarly strong linear relationship has been found between the exit velocity of the exhaust gas and the downward dispersion (Fig. 7 and Table 2) with regression R² of 1.00 for cases with and without obstacles. It is, however, a negative dependence, because higher exit velocities transport the plume into higher altitudes, and consequently the downward dispersion is lower. The effective range is much smaller than for the wind speed with only 3.7 % with and 2.1 % without obstacles, respectively.

Figure 7 Dependence of the downward dispersion D on different exit velocities v_exit with and without the obstacle effect.

3.1.2 Effect of wind direction

The strength of the downward dispersion was investigated depending on different wind directions in relation to the orientation of the ship. Frontal wind (angle of 0^∘) hits the short side of the vessel, which has a width of 30 m, whereas lateral wind (angle of 90^∘) has to be lifted over the 246 m length of the ship. Therefore, a stronger distortion of the flow during lateral wind has been observed.

The downward dispersion correlates linearly with the cosine of the flow angle ϕ (Fig. 8). A regression coefficient R² of 0.98 was calculated. At default settings, downward dispersion ratios of 7.0 % and 16.6 % was found under frontal and lateral wind conditions, respectively. This results in an effective range of 9.6 %. The corresponding downward dispersion under a no-obstacle condition is 2.3 %. There is no effective range for no-obstacle conditions, because here a single symmetrical stack is assumed, where the downward dispersion values are the same for both frontal and lateral winds. However, very small differences between these conditions can occur during the modeling (see Table C1), which result from asymmetry in the numerical grid.

Figure 8 Dependence of the downward dispersion D on the cosine of different wind flow angles (ϕ) towards the ship. An angle of 0^∘ indicates the frontal wind; 90^∘ indicates the lateral wind.

3.1.3 Effect of exhaust plume temperature

The exhaust plume temperature depends on technical parameters like the engine power and the use of a heat exchanger, and therefore a range of possible temperatures (200–400 ^∘C) was investigated. Figure 9 presents results of the single linear regression for the downward dispersion at varying exhaust temperatures with and without the obstacle effect.

Figure 9 Dependence of the downward dispersion D on different exhaust temperatures T_exh with and without the obstacle effect.

Once again, a strong linear relationship with correlation coefficients R² of 0.98 and 0.99 was found for results with and without the ship, respectively. At higher exhaust temperatures, the plume reaches higher altitudes by convective upward movement, which results in lower downward dispersion ratios. The effective range under default settings is 6.9 % with and 2.9 % without the obstacle effect (Table 2).

3.1.4 Effect of atmospheric stability

The effect of atmospheric stability Γ on the downward dispersion was investigated in a range from unstable (−1.2 K × 100 m⁻¹) to very stable (+0.5 K × 100 m⁻¹) vertical temperature gradients. Under default settings, linear regression resulted in correlation coefficients of R²=0.90 and 0.94 with and without the ship, respectively (Fig. 10a). Since the R² coefficient was low compared to the other investigated input parameters, a linear dependence would deliver poorer results for this parameter. Therefore, a quadratic dependence was calculated as well.

Figure 10 Dependence of the downward dispersion D on different vertical temperature gradients Γ with and without the obstacle effect. Linear regressions of the downward dispersion against Γ (a) and sgn(Γ)Γ² (b) are shown.

Since the square of a negative vertical temperature gradient would result in a positive value, a sign function was applied. The mathematical expression is

$$\begin{array}{}\text{(3)}& \mathrm{sgn}\left(x\right):=\left\{\begin{array}{ll}-\mathrm{1}& \text{if}x<\mathrm{0}\\ \mathrm{0}& \text{if}x=\mathrm{0}\\ \mathrm{1}& \text{if}x>\mathrm{0}\end{array}\right.\end{array}$$

Then, the correlation between downward dispersion and sgn(Γ)Γ² is calculated (Fig. 10b). It shows better agreement in the cases considering obstacle effects (R² = 0.99) and slightly better agreement in cases without the ship (R² = 0.96), as well. It is a negative correlation because higher temperature gradients correspond to a higher stability, which thermodynamically prevents the plume from dispersing vertically and therefore lowers the downward dispersion ratio. The effective range of the temperature gradient on the downward dispersion is 6.6 % for ship cases and 3.8 % for stack-only cases.

3.2 Result of the multiple regression

Multiple regression is performed according to the equations in Appendix C2. The downward dispersion ratio depends linearly on all investigated input parameters, their cosine (in the case of the angle of wind direction) or their squares (in the case of atmospheric stability). With that in mind, a training dataset for the multiple regression was created. Here, all independent input parameters are varied at the same time (but in the given range), and the downward dispersion ratio is calculated with MITRAS. For a set of 39 different combinations (Table C1) of input parameters with the obstacle effect and 27 without, the estimation coefficients ${\widehat{\mathit{\beta }}}_i$ for individual parameters i (wind speed, exit velocity, etc.) are calculated with the multiple regression. The number of simulated cases without obstacle effects is lower, because in these cases the wind direction has been varied which will not show differences in the case of stack-only conditions. The resulting formulae for the parameterization read

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& D\left[\mathit{\%}\right]=\mathrm{13.03}+\mathrm{3.45}{v}_{\mathrm{wind}}-\mathrm{1.01}{v}_{\mathrm{exit}}-\mathrm{0.026}{T}_{\mathrm{exh}}\\ & -\mathrm{3.81}\mathrm{sgn}\left(\mathrm{\Gamma }\right){\mathrm{\Gamma }}^{\mathrm{2}}+\mathrm{6.13}\cos \left(\mathit{\varphi }\right),\end{array}\end{array}$$

with the ship and

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& D\left[\mathit{\%}\right]=\mathrm{4.55}+\mathrm{178}{v}_{\mathrm{wind}}-\mathrm{64}{v}_{\mathrm{exit}}-\mathrm{0.018}{T}_{\mathrm{exh}}\\ & -\mathrm{3.40}\mathrm{sgn}\left(\mathrm{\Gamma }\right){\mathrm{\Gamma }}^{\mathrm{2}},\end{array}\end{array}$$

without the ship (i.e., stack-only).

Here, v_wind and v_exh are given in [m s⁻¹], T_exh in [^∘C], Γ in [K/100 m] and ϕ in [^∘], where 0^∘ refers to frontal wind and 90^∘ to lateral wind.

3.3 Bootstrapping

A bootstrapping procedure is applied to estimate how well the parameterization can represent the model data. For this purpose, downward dispersion ratios were calculated with the parameterization formulae (Eqs. 4 and 5) and compared to the original MITRAS results for all investigated cases and ranges. The individual parameterization results are listed in Table C1. Table 3 gives the overall results of the bootstrapping procedure.

Table 3 Results of the bootstrapping procedure for cases with and without considering the ship-induced obstacle effect.

With a mean absolute error of 1.9 ± 1.6 % for cases with the ship and 1.2 % ± 0.9 % without the ship, the parameterization delivers very similar results to the model runs. The maximum absolute errors were found to be 6.1 % in cases with the ship and 4.0 % in cases without the ship.

3.4 Assessment of the obstacle effect

Another aim was to investigate under which conditions the strongest downward dispersion occurs and which effect the consideration of the obstacle has on the downward dispersion.

From the single-parameter regressions, it is assumed that the strongest downward dispersion occurs at high wind speed (15 m s⁻¹) with lateral wind (90^∘), low exit velocity (4 m s⁻¹), low plume temperature (200 ^∘C) and during unstable atmospheric conditions (−1.2 K × 100 m⁻¹). This is displayed in Fig. 11 with the ship as an obstacle (panel a) and under stack-only conditions (panel b).

Figure 11 Visualization of the obstacle effect in MITRAS. Examples for lateral wind with v_wind=15 m s⁻¹, v_exit= 4 m s⁻¹, T_exh=200 ^∘C and $\mathrm{\Gamma }=-\mathrm{1.2}$ K × 100 m⁻¹ are presented. Concentration fields are displayed with the ship pointing towards the viewer (a) and under stack-only conditions (b).

The calculated downward dispersion ratio for this condition is 54.9 % and 31.1 % with and without the obstacle effect, respectively. This means that a significant proportion of nearly 25 % of the emission can be dispersed downwards only by taking into account the turbulence caused by the ship.

3.5 Discussion of limitations

Despite efforts to represent real conditions as best as possible, the results are subject to a few limitations or uncertainties that will be discussed in the following section. For a general comparison of our MITRAS results with a common dispersion model, see Appendix D.

One factor that is not considered in this study is relative humidity. Here, a distinction must be made between the relative humidity of the ambient air and the relative humidity of the exhaust. By using the Lagrangian concept based on the so-called projected area entrainment (Lee and Cheung, 1990), Affad et al. (2006) stated that the relative humidity of the ambient air has only a slight impact on the plume rise, diameter and temperature for values between 20 % and 90 %. It can have an impact on particle growth, but as this study focuses on a passive gaseous tracer, this effect is neglected. On the other hand, the humidity of the exhaust gas might have a larger impact on the plume rise. Since water vapor has a lower density than air, an exhaust gas mixture of high humidity will show a stronger plume rise. Furthermore, as the gas will quickly condense, it will release latent heat and rise further. However, the database on humidity of ship exhaust is sparse. It could play a role in the case of vessels using a scrubber to wash out SO₂ from the exhaust. During this process, the exhaust is cooled down significantly and therefore will show a weaker plume rise (Murphy et al., 2009). It is unclear if the additional buoyancy can compensate for the lower exhaust temperatures. Due to these uncertainties and lack of data, the relative humidity has not been included in this study.

Second, the emission is assumed to occur in the grid cell above the stack, which has a size of 2 m × 2 m × 2 m. This corresponds to a stack with a square cross section of 4 m² and is a limitation connected to the chosen grid size. Real stacks are usually round and have a smaller diameter. However, Bai et al. (2020) reported about plume modeling for container vessels with measured funnel diameters in the range of 1.38 to 3.0 m, so the selected value is inside a reasonable range. Furthermore, many ships operate multiple smaller stacks that might in sum lead to a similar exhaust behavior. Based on the selected diameter, the real exit velocity could differ slightly. However, by comparing the effective ranges for exit velocity against all other input factors (Table 2), it can be seen that this parameter has the smallest overall impact of the downward dispersion, and therefore this uncertainty factor has a low impact on the overall performance.

Another assumption was that the ship has been considered as a non-moving source, i.e., a hoteling ship. However, the results can be applied to a moving ship by calculating the vector sum of the wind and the vessel speed. It is difficult to account for complex maneuvers, though, as the resulting wind vector may change quickly and the technical conditions like exhaust temperature and exit velocity may also vary with the speed of the ship.

The shape of the vessel and the location of the stack are additional parameters that can influence the exact value of downward dispersion. Parameterizing them is beyond the scope of this study as the shape was chosen to investigate the average effect of a cruise-ship-sized vessel on pollutant concentrations close to ground. However, to get an impression, an exemplary comparison of MITRAS results for a cruise ship and container vessel can be found in Table S2 in the Supplement.

The chosen model surface is water but assuming a hoteling ship, the land surface effects may play a role for the dispersion. This effect has not been part of this study, as this is a highly variable parameter that depends on the structure of the harbor, the city and the orography. These effects need to be covered by a larger-scale model.

4 Conclusions

A ship plume modeling study was performed with the microscale numerical model (MITRAS) to investigate the downward dispersion of the exhaust in close proximity to a modeled cruise ship (i.e., in the momentum-driven regime). A set of 39 different scenarios with varying meteorological and technical input parameters was analyzed. A multiple regression algorithm was used to estimate a parameterization function for the downward dispersion. This parameterization has been tested against the MITRAS model results through a bootstrapping procedure.

From single-parameter regressions, a positive linear relationship of the downward dispersion from wind speed and negative linear relationships from exit velocity, plume temperature and the cosine of the angle of wind direction were found. The downward dispersion ratio was larger in the case of lateral wind than in the case of frontal wind. In the case of atmospheric stability, the downward dispersion showed a squared dependence from the vertical temperature gradient multiplied by the sign function. From all these input parameters, the wind speed shows the largest effect on the downward dispersion in the investigated range (2–15 m s⁻¹).

A comparison of the model results and the parameterization from multiple regression shows a good agreement with a mean absolute error of 1.9 ± 1.6 % for cases with the ship and 1.2 ± 0.9 % without the ship. For the case of strongest downward dispersion, the difference was calculated between downward dispersion with (54.9 %) and without considering the obstacle effect (31.1 %), which was almost 25 %. This shows how important it is to consider the effects of the downward dispersion in the momentum-driven regime when one wants to evaluate the air pollution situation in harbor areas.

The parameterization functions can also be used for container ships of a similar size. It may also be applied to different emission situations like industrial stacks.

In a future study, other plume parameters will be derived from the vertical concentration profiles in a similar way to the downward dispersion. This includes the height of the plume axis and the shape of the vertical plume profile, which may deviate from the often-assumed Gaussian distribution. These results can further be used in a city-scale model, which only calculates the plume dispersion inside the buoyancy-driven regime.

Appendix A: Wind speed in Hamburg

Input data for wind speed data were used from the Hamburg weather mast, provided by the Integrated Climate Data Center (ICDC) (ICDC, 2020). The weather mast is positioned at a meteorological measurement station in Billwerder, Hamburg ($\mathrm{53}{}^{\circ }{\mathrm{31}}^{\prime }{\mathrm{09.0}}^{\prime \prime }$ N, $\mathrm{10}{}^{\circ }{\mathrm{06}}^{\prime }{\mathrm{10.3}}^{\prime \prime }$ E). Hourly data from all of 2018 were statistically analyzed at five different measurement heights (Fig. A1, Table A1).

Figure A1 Boxplots of hourly wind speed data for the year 2018 from the Hamburg weather mast. Red lines indicate median values; lower and upper whiskers end at 5th and 95th percentile, respectively.

Table A1 Statistical data on hourly wind speed values [m s⁻¹] from the Hamburg weather mast for the year 2018.

Appendix B: Plume temperature

Figure B1 presents results for maximum temperatures in the MITRAS domain for one case with the highest temperature (400 ^∘C). Ambient temperatures (15 ^∘C) are reached at a horizontal distance of approximately 100 m from the stack.

Figure B1 Results for calculated maximum temperatures in the MITRAS domain in distance downwind from the stack for a case of a 400 ^∘C plume with an ambient temperature of 15 ^∘C. The values closer to the stack are not shown for clarity. The thermic effect on the plume rise due to the hot exhaust gas vanishes at around 100 m distance from the stack.

Appendix C: Regressions

This section describes the general application of linear and multiple regression on the model results.

C1 Single linear regression

A simple approach to estimate a target variable Y (e.g., the downward dispersion) from a single independent variable X (e.g., the wind speed) is a linear regression in the form of

$$\begin{array}{}\text{(C1)}& \widehat{Y}=\widehat{{\mathit{\beta }}_{\mathrm{0}}}+\widehat{{\mathit{\beta }}_{\mathrm{1}}}X,\end{array}$$

where β₀ and β₁ are the ordinate axis intersection and the slope, respectively, and the circumflex (∧) describes an estimated parameter. $\widehat{{\mathit{\beta }}_{\mathrm{0}}}$ and $\widehat{{\mathit{\beta }}_{\mathrm{1}}}$ are calculated with the least squares method, minimizing the quadratic deviation between model result values Y_i and estimated values $\widehat{Y_i}$. The required function Q reads

$$\begin{array}{}\text{(C2)}& \begin{array}{rl}& Q\left(\widehat{{\mathit{\beta }}_{\mathrm{0}}},\widehat{{\mathit{\beta }}_{\mathrm{1}}}\right)={\sum }_{i=\mathrm{1}}^n{\left(Y_i-\widehat{Y_i}\right)}^{\mathrm{2}}=\\ & {\sum }_{i=\mathrm{1}}^n{\left(Y_i-\widehat{{\mathit{\beta }}_{\mathrm{0}}}-\widehat{{\mathit{\beta }}_{\mathrm{1}}}{X}_i\right)}^{\mathrm{2}}.\end{array}\end{array}$$

Minimizing is done by applying partial derivations from Q to $\widehat{{\mathit{\beta }}_{\mathrm{0}}}$ and $\widehat{{\mathit{\beta }}_{\mathrm{1}}}$. This results in

$$\begin{array}{}\text{(C3)}& {\displaystyle }\widehat{{\mathit{\beta }}_{\mathrm{0}}}=\overline{Y}-\widehat{{\mathit{\beta }}_{\mathrm{1}}}\overline{X},\text{(C4)}& {\displaystyle }\widehat{{\mathit{\beta }}_{\mathrm{1}}}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n\left(\left(X_i-\overline{X}\right)\left(Y_i-\overline{Y}\right)\right)}{{\sum }_{i=\mathrm{1}}^n{\left(X_i-\overline{X}\right)}^{\mathrm{2}}}},\end{array}$$

where $\overline{X}$ and $\overline{Y}$ are the mean values of the corresponding dataset.

C2 Multiple regression

The variable Y can depend on more than one independent input variable (X₁, X₂, …, X_p). Then, a multiple regression can be applied, and in the case of linear dependencies, the corresponding regression is called multilinear. The multilinear estimation for $\widehat{Y}$ reads

$$\begin{array}{}\text{(C5)}& \widehat{Y}=\widehat{{\mathit{\beta }}_{\mathrm{0}}}+\widehat{{\mathit{\beta }}_{\mathrm{1}}}{X}_{\mathrm{1}}+\widehat{{\mathit{\beta }}_{\mathrm{2}}}{X}_{\mathrm{2}}+\mathrm{\dots }+\widehat{{\mathit{\beta }}_p}{X}_p.\end{array}$$

Again, the minimum distance between Y_i and $\widehat{Y_i}$ can be calculated by the least squares method, similar to the case of linear regression, by minimizing the function Q:

$$\begin{array}{}\text{(C6)}& \begin{array}{rl}& Q\left(\widehat{{\mathit{\beta }}_{\mathrm{0}}},\widehat{{\mathit{\beta }}_{\mathrm{1}}},\widehat{{\mathit{\beta }}_{\mathrm{2}}},\mathrm{\dots }\widehat{{\mathit{\beta }}_p}\right)={\sum }_{i=\mathrm{1}}^n\\ & {\left(Y_i-\widehat{{\mathit{\beta }}_{\mathrm{0}}}-\widehat{{\mathit{\beta }}_{\mathrm{1}}}{X}_{i,\mathrm{1}}-\widehat{{\mathit{\beta }}_{\mathrm{2}}}{X}_{i,\mathrm{2}}-\mathrm{\dots }-\widehat{{\mathit{\beta }}_p}{X}_{i,p}\right)}^{\mathrm{2}}.\end{array}\end{array}$$

However, as this can lead to complicated expressions of $\widehat{\mathit{\beta }}$, one can make use of a matrix representation.

$$\begin{array}{}\text{(C7)}& \begin{array}{rl}& \left(\begin{array}{c}\begin{array}{c}\widehat{Y_{\mathrm{1}}}\\ \widehat{Y_{\mathrm{2}}}\end{array}\\ \begin{array}{c}\mathrm{⋮}\\ \widehat{Y_n}\end{array}\end{array}\right)=\left(\begin{array}{c}\begin{array}{c}\begin{array}{ccc}\mathrm{1}& X_{\mathrm{11}}& \begin{array}{ccc}{X}_{\mathrm{12}}& \mathrm{\cdots }& X_{\mathrm{1}p}\end{array}\end{array}\\ \begin{array}{ccc}\mathrm{1}& X_{\mathrm{21}}& \begin{array}{ccc}{X}_{\mathrm{22}}& \mathrm{\cdots }& X_{\mathrm{2}p}\end{array}\end{array}\end{array}\\ \begin{array}{c}\begin{array}{ccc}\mathrm{⋮}& \mathrm{⋮}& \begin{array}{ccc}\mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\end{array}\end{array}\\ \begin{array}{ccc}\mathrm{1}& X_{n\mathrm{1}}& \begin{array}{ccc}{X}_{n\mathrm{2}}& \mathrm{\cdots }& X_{np}\end{array}\end{array}\end{array}\end{array}\right)\\ & \cdot \left(\begin{array}{c}\begin{array}{c}\widehat{{\mathit{\beta }}_{\mathrm{0}}}\\ \widehat{{\mathit{\beta }}_{\mathrm{1}}}\end{array}\\ \begin{array}{c}\mathrm{⋮}\\ \widehat{{\mathit{\beta }}_p}\end{array}\end{array}\right).\end{array}\end{array}$$

By using both the transpose (^T) and the invert (⁻¹) operators, the equations can be transformed to a general solution for $\widehat{\mathit{\beta }}$:

$$\begin{array}{}\text{(C8)}& \left(\begin{array}{c}\begin{array}{c}\widehat{{\mathit{\beta }}_{\mathrm{0}}}\\ \widehat{{\mathit{\beta }}_{\mathrm{1}}}\end{array}\\ \begin{array}{c}\mathrm{⋮}\\ \widehat{{\mathit{\beta }}_p}\end{array}\end{array}\right)={\left({\mathbf{X}}^T\mathbf{X}\right)}^{-\mathrm{1}}\cdot {\mathbf{X}}^T\cdot \left(\begin{array}{c}\begin{array}{c}\widehat{Y_{\mathrm{1}}}\\ \widehat{Y_{\mathrm{2}}}\end{array}\\ \begin{array}{c}\mathrm{⋮}\\ \widehat{Y_n}\end{array}\end{array}\right).\end{array}$$

Table C1 presents an overview of the results of the multiple regression.

Table C1 Data table for regression analyses. v_{wind, in} refers to the input wind speed at the top model layer; v_{wind, stack} refers to the wind speed at stack height. Further input data are exit velocity (v_exit), exhaust temperature (T_exh) and wind direction (ϕ), with 0^∘ referring to frontal and 90^∘ to lateral wind. Results are given for downward dispersion with and without the obstacle effect (D and D_stack-only). D_par refers to results of the parameterization. The bold values in line no. 8 correspond to the default settings. Values in brackets in D_{par, stack-only} were not included in the multiple regression, because in these cases only the wind direction was changed, which does not affect stack-only results.

Appendix D: Comparative assessment of the model results

To get an impression of similarities and differences of MITRAS results with a common dispersion model, the results are compared under similar conditions with results from the integral plume model IBJpluris (Janicke and Janicke, 2001), which can be used to describe the plume dispersion in the momentum-driven regime. IBJpluris calculates average plume properties like concentration and temperature along the plume centerline and applies a circular Gaussian dispersion around this central axis. IBJpluris does not account for obstacle-induced turbulence effects and is therefore only compared to stack-only conditions in MITRAS.

Since the primary output of IBJpluris is the plume centerline and not the downward dispersion, a similar centerline height for MITRAS was calculated to compare the plume behavior. Therefore, the centerline in MITRAS h_{center, MITRAS} is defined as the median height of the plume mass (i.e., 50 % of the plume mass lies below and 50 % lies above). It is calculated at the same distance as downward dispersion for a column of 100 m × 100 m (see Fig. 4). Since this is an average of values between a distance of 100 to 200 m, IBJpluris centerline heights were calculated at a distance of 100 to 200 m as well (h_{center, Pluris}). Table S1 gives an overview of the comparison. Δh_MITRAS and Δh_Pluris are the differences between plume height (52 m) and centerline height for MITRAS and IBJpluris calculations, respectively. Their minimum difference is given by $min\left(\right|\mathrm{\Delta }{h}_{\text{MITRAS}}-\mathrm{\Delta }{h}_{\text{Pluris}}\left|\right)$ in Table S1, which represents the closest similarity of both models. Results are given at default settings and selected conditions to compare effects of input parameters. For all selected cases, MITRAS calculates larger centerline height values than IBJpluris. The lowest differences occur at low wind speed, low exhaust temperature and very stable conditions. The strongest differences of over 20 m occur for cases of low exit velocity and high exhaust temperature.

By calculating effective ranges, $r_i=\mathrm{\Delta }|min(|\mathrm{\Delta }{h}_{\text{MITRAS}}-\mathrm{\Delta }{h}_{\text{Pluris}}|{)}_{i,max}-min\left(\right|\mathrm{\Delta }{h}_{\text{MITRAS}}-\mathrm{\Delta }{h}_{\text{Pluris}}\left|{)}_{i,min}\right|$ for a certain input parameter i, one can evaluate which input parameter causes the highest discrepancy between the models. For example, changing the wind speed only results in an effective range of 1 m, while temperature and stability changes both show effective ranges of 10 m.

The higher plume rise in MITRAS is consistent with the interaction of the hot plume with the ambient air. MITRAS accounts for the change in the thermodynamic field, and the heat balance equation creates an additional buoyancy which is not considered in simpler approaches. This explains the high effective range for temperature and stability changes.

This shows that the results for stack-only conditions are reasonable and that MITRAS provides a more complex improvement over simple Gaussian approaches in the near field.