Diagnosing topographic forcing in an atmospheric dataset: The case of the North American Cordillera

It is well known from modelling studies that surface topography influences the large‐scale atmospheric circulation and that several model biases are associated with incorrect representation of topography. The textbook explanation of topographic effects on large‐scale circulation appeals to the theoretical relationship between surface forcing and vortex stretching along trajectories in single‐layer models. The goal of this study is to design and use a simple diagnostic of the large‐scale forcing on the atmosphere when air is passing over topography, directly from atmospheric fields, based on this theoretical relationship. The study examines the interaction of the atmosphere with the North American Cordillera and samples the flow by means of trajectories during Northern Hemisphere winter. We detect a signal of topographic forcing in the atmospheric dataset, which, although much less distinct than in the theoretical relationship, nevertheless exhibits a number of expected properties. Namely, the signal increases with latitude, is usually stronger upslope than downslope, and is enhanced if the flow is more orthogonal to the mountain ridge, for example during periods of positive Pacific–North American index (PNA). Furthermore, a connection is found between an enhanced signal of topographic forcing downslope of the North American Cordillera and periods of more frequent downstream European blocking.


INTRODUCTION
The zonal asymmetry of the climatological atmospheric circulation on Earth is primarily induced by zonal asymmetries of the lower boundary conditions.The impact of topography, land-sea contrasts and sea-surface temperature (SST) patterns on the circulation is evident both locally and remotely (Charney and Eliassen, 1949;Inatsu et al., 2002).Topography induces flow anomalies because air is forced either to ascend or to flow around an obstacle, depending on atmospheric stability and the large-scale structure of the topography (Valdes and Hoskins, 1991;Lott and Miller, 1997).In the simplest and traditional way, the topographic impact can be understood by solving the vorticity equation in a Lagrangian framework and focusing on the behaviour of the forcing by the stretching/divergence term (e.g.Holton, 2004).In the presence of a steady flow, topography initiates stationary Rossby waves (Charney and Eliassen, 1949), which communicate the impact of the topography downstream, both horizontally and vertically (e.g.Hoskins and Karoly, 1981;Held, 1983;Held et al., 2002).
In atmospheric models, both resolved and unresolved topographic processes need to be captured to ensure correct representation of the complete effect of topography on the system.The impact of topography is partly included in physical parametrizations that modify the circulation via low-level drag and gravity-wave drag (e.g.Palmer et al., 1986;Lott and Miller, 1997;Sandu et al., 2016).It is known that the representation of topography is crucial for model fidelity of phenomena downstream, for example the North American Cordillera and specifically the Rocky Mountains for the Atlantic storm track (Brayshaw et al., 2009;Pithan et al., 2016).The representation of topography, and potentially of the Rocky Mountains, also modulates atmospheric blocking over the Atlantic and Europe (Jung et al., 2012;Berckmans et al., 2013;Pithan et al., 2016).
From the perspective of model biases it is desirable to detect if topographic forcing, both resolved and parametrized, is represented correctly.The importance of topography for flow phenomena and related biases is typically analysed by "mountain -no mountain" experiments (e.g.Held et al., 2002;White et al., 2017) or by comparing atmospheric circulation with and without topographic parametrizations (e.g.Pithan et al., 2016).However, in such experiments the subsequent change of the circulation induces adjustments, for example in diabatic heating and nonlinear interactions, implying that a combination of several factors is responsible for the total alterations in the patterns of the atmospheric circulation.As a result, the footprint of the direct topographic forcing might be masked.To our knowledge, no study has previously tried to detect a signal of the topographic forcing expected from theory directly from the atmospheric fields themselves.
In this work we evaluate if and to what degree a signal of the topographic forcing over large-scale elevated areas, as suggested by the vorticity tendency equation through the stretching term, can be found in an atmospheric dataset.Furthermore, we investigate if the atmospheric flow takes a path that is more often over or around the highest elevation of the topography.The former is mostly assumed for latitudinally broad mountains (Valdes and Hoskins, 1991).
In line with the traditional picture of following a column of air while it is interacting with the topography (e.g.Holton, 2004), this study applies Lagrangian sampling before investigating the impact of the topography on the large-scale circulation.The dataset used consists of about 500 trajectory sets (around 2 ⋅ 10 5 individual trajectories) calculated from wind fields of the ERA-Interim reanalysis (Dee et al., 2011), which we use to characterise the flow over topography in the real atmosphere.For the subsequent calculation of the topographic forcing the Lagrangian samples are aggregated together.To illustrate the method, we focus on the atmospheric interactions with the extended Rocky Mountains (North American Cordillera) in Northern Hemisphere midlatitudes.The diagnostic tool developed here highlights how different flow regimes are linked to the strength of the topographic forcing signal.Section 2.1 develops the diagnostic, Sections 2.2 and 2.3 discuss the other aspects of the methodology, Section 3 presents the results, and Section 4 summarises the results and discusses some of the limitations of the analysis.

Topographic forcing in shallow-water dynamics
To motivate our diagnostic, we start by recalling the impact of localised topography on an impinging zonal flow within the shallow-water equations.Figure 1 depicts the quasi-steady-state solution for a particular set of parameters, as a reminder of the qualitative picture (a detailed description of the model set-up can be found in Appendix A).A trajectory (panel a) exhibits the characteristic pressure ridge over the topography and the downstream lee trough (Figure 1a,b), in line with dynamic theory for large-scale, quasi-steady flows with Froude number on the order of the Rossby number and thus consistent with quasi-geostrophic scaling.The horizontal motion shown can be explained by the shallow-water vorticity equation which, in pressure coordinates and in a Lagrangian frame of reference with the Lagrangian derivative D/Dt, is given by: We aim to apply this theory to large-scale motions and therefore assume quasi-geostrophic scaling.The stretching term (term on the right-hand side) is the product of the horizontal divergence  and the absolute vorticity, composed of the relative vorticity  and the Coriolis term f .F I G U R E 1 Characteristics of the quasi-steady-state solution of a shallow-water model describing flow over a Gaussian ridge (black contour lines, at 100 m and 1,000 m).The panels show the (a) depth of flow, (b) relative vorticity and (c) horizontal divergence, in a -plane channel, in coloured contours.Contour lines of positive (negative) values are depicted in full (dashed) lines and the zero line is omitted.In panel (a) a trajectory (streamline in quasi-steady-state), which is initialised at x = 0 km and y = 3,500 km, is visualised as an orange line.In the y-direction the central 4,000 km of the channel are shown and all fields are smoothed by applying a 1-2-1 filter three times to balance weak numerical dissipation It captures the compression of the flow upslope (positive horizontal divergence) and stretching downslope of the topography (negative horizontal divergence, Figure 1c).The upslope compression translates to a negative relative vorticity tendency because the absolute vorticity,  + f , has the sign of f ; this is true in the model, and almost everywhere in the real atmosphere.In the same way, downslope stretching imparts a positive vorticity tendency.The interaction with the topography also induces a stationary Rossby wave downstream of the ridge.
The topographic forcing is approximately communicated by a linear relationship between the surface topography h and the divergence .In the shallow-water case, this relationship can be computed directly using Equation 1 and the height equation of the shallow-water model which follows from mass conservation.In the following, H 0 is the depth of the flow upstream of the topography and h top the component of the flow depth along streamlines.On an f -plane, the upstream value f /H 0 is uniform in Equation 3 and it can be used to substitute for H on the right-hand side of Equation 2. On a -plane, the variation of the upstream value of f /H 0 sampled by a typical latitudinal displacement of the trajectory can be neglected for mountains with characteristic cross-mountain scale much smaller than √ U  ≈ 10 6 m.This condition is not necessarily strictly satisfied for the mountain chains of interest, but for a potential diagnostic relation we will make the approximation and perform the substitution.Then the steady state of Equation 2 can be written in the approximate form: Using the geostrophic vorticity  = g f ∇ 2 h top in the PV-equation (Equation 3) gives an elliptic equation for h top , namely . From this it can be derived that for horizontal scales much smaller than √ gH 0 f ≈ 3,000 km the order of magnitude of h top is much less than that of h, so that the first term on the left-hand side of Equation 4 can be neglected compared with the second: (5) This is the approximate diagnostic relationship between the surface forcing f sfc = ) and the stretching term that forms the basis of the diagnostic developed in this article.It should be noted that a simple damping term − r can be added to the vorticity equation (Equation 1), following Held (1983), and could be carried through to this diagnostic relation.Equation 5 is advantageous compared to Equation 4 because we seek a simple, robust empirical measure of the signal of the surface topographic forcing that can be evaluated in a complex atmospheric dataset, where the component of the flow interacting directly with the topography is in general unknown.The description of the surface forcing is also in line with the linearised boundary condition for a broad mountain w ′ = u h ′ x , where the prime indicates deviations from the zonal mean (Valdes and Hoskins, 1991).
Figure 2 shows the idealised results in the phase space defined by the stretching term and the surface forcing f sfc for two experiments.Experiment (a) is the default set-up, for which results are shown in Figure 1.A second experiment shows the impact of increased dissipation through the use of a longer time step (Δt = 60 s instead of Δt = 3 s), while keeping the shape of the topography the same.The second experiment thus slightly deviates from the relationship described in Equation 5.It presents an example of how additional processes affect the topographic forcing, as measured by the relation between surface forcing and stretching term.
In the phase space, the upslope evolution along the trajectory corresponds to the branch extending to the bottom-right quadrant.Passing the peak of the idealised mountain the surface forcing becomes negative and the evolution continues in the top-left quadrant of the phase space.A slope of around −10 −8 m −1 ⋅s −1 in both set-ups is in line with Equation 5, based on f = 10 −4 s −1 and H 0 = 10 4 m.
In the second experiment (panel b) a broadening of the curve can be explained by the fact that the longer time step increases the implicit numerical dissipation, which introduces a phase shift in the divergence field relative to the topography by horizontally smoothing out the H-anomalies.The inverse of the effective dissipation time-scale can be calculated for this experiment from the ratio of the stretching term and relative vorticity at the top of the ridge (Held, 1983), based on Equation 1 in steady-state conditions with an additional damping term −r on the right-hand side, as r = − 1b and 2b).Thus, the effective dissipation time-scale is on the stronger side compared to estimates by Charney and Eliassen (1949) of r −1 ≈ 6 days and r −1 ≈ 12 days for barotropic and baroclinic currents, respectively.Nevertheless, the linear relationship is essentially preserved.We cannot expect the real atmosphere to exhibit such a clean linear relation as in the idealised model, not least because the effective flow depth H 0 , which affects the slope of the relation, will depend on the flow conditions.Thus we seek a diagnostic that is sensitive to the qualitative behaviour predicted by the theory, but insensitive to quantitative details.The theoretical relationship predicts a negative correlation between surface forcing and F I G U R E 2 Results from the shallow-water experiments shown in the phase space defined by stretching term and surface forcing f sfc .Small grey points represent results from all grid points of the steady-state solution and larger coloured points highlight points along the central trajectory (as in Figure 1a).The shading of the coloured points varies with the evolution of the trajectory as darker points correspond to locations further east.The two experiments are described in Section 2.1.In panel (a) the q i -labels map the number of trajectory points q i to their respective quadrants, which is relevant for Equation 6stretching, namely an increased population of points in the bottom-right and top-left quadrants of the phase space as compared to the other two quadrants.Comparison of the two experiments in Figure 2 shows that in both a negative correlation is present, although slightly modified in panel b through the additional damping term.Thus, we can evaluate the predicted relation categorically, and qualitatively capture a signal of the topographic forcing by determining the asymmetry of the number of points q i in quadrants i ∈{1, 2, 3, 4} in phase space, since a random distribution would be symmetric.The mapping of the total number of trajectory points q i to quadrants i is shown in Figure 2a.The effectiveness of a predicted categorical relationship is measured by the odds ratio, which in the the horizontal wind at 700 hPa during one period of positive PNA index.Longer and darker arrows indicate higher wind speeds.Trajectories are released from the bold box at 140 • W, indicating the initial curtain of trajectory points.Trajectories are only analysed between the two bold-dashed meridional lines.Grey contour lines (every 500 m) represent the North American Cordillera and become darker with increasing elevation upstream and downstream regions can be calculated as respectively.Values larger than unity we take as evidence of positive topographic forcing.

Trajectory tool and topographic forcing analysis
The diagnostic must be applied only where air parcels are interacting with topography.We determine this through trajectories.The Reading Offline Trajectory Model (ROTRAJ) uses a 4th-order Runge-Kutta scheme to calculate trajectories based on the three-dimensional gridded wind field (Methven et al., 2001).Information on the wind and other variables at each trajectory point is obtained by interpolating bi-linearly in time and horizontal space, with cubic Lagrange interpolation in the vertical.To prevent trajectories from crossing the surface, calculations are done on vertical -levels, which are the hybrid-pressure terrain-following coordinate of the Integrated Forecasting System (IFS) model (Simmons and Burridge, 1981).For increased accuracy the native IFS model output of spectral fields on -levels is used as input to the trajectory model.The tool allows calculation of forward trajectories and back trajectories.
A set of trajectories is initialised along 140 • W at 61 latitude points between 30 and 60 • N (Figure 3).Trajectories are released at 10 pressure levels of which six are used for the default analysis (850 to 600 hPa, spaced by 50 hPa).From the initial curtain of trajectory points (bold box in Figure 3), forward trajectories are calculated over 4 days, using 6-hourly ERA-Interim reanalysis data (Dee et al., 2011) with an integration step of 1 h.This set-up is repeated consecutively every 4 days to obtain 546 winter (December to February) trajectory sets for the period 1980-2005 where each set contains one 4-day analysis.This means that in total around 2 ⋅ 10 5 trajectories are analysed.As the focus of this study is on the interaction of the flow with the North American Cordillera, analysis is limited to points between 90 and 140 • W (thick meridional lines in Figure 3).In addition, trajectory points above 300 hPa are discarded to exclude events that are likely not influenced by the topography.Several parameters are stored along the trajectories, namely the horizontal position, pressure, horizontal wind, Ertel PV, divergence and absolute vorticity.The sensitivity of the topographic forcing signal to some of the data selection choices discussed here is given in Appendix B.
Zonal and meridional gradients of surface topography are calculated based on the topography adopted in ERA-Interim.The surface forcing term f sfc is calculated from the gradient of the topography and is then interpolated to the individual trajectory points.It should be noted that the horizontal wind at the height of each trajectory point is used instead of the wind at a constant height or close to the surface.This approach is chosen to simplify the analysis and assumes implicitly that the atmosphere is barotropic, which is not in general true.
Once the set of trajectory points is compiled, the stretching term in pressure coordinates (calculated from native -coordinates as described in Appendix C) can be related to the surface forcing for each point to determine the topographic forcing .In the following, calculations of  (Equation 6) are either based on all points in the phase space or the data are limited to only include signals stronger than a certain threshold.For the latter, points are discarded if their stretching term is smaller than one standard deviation of the distribution of all trajectory sets and their surface forcing smaller than half of the standard deviation.This is done as a simple way of setting a noise threshold, since data errors would preferentially affect small values of the terms.The thresholds do not bias the diagnostic since they are symmetric in the phase space.The second approach of restricting the analysis to points with stronger signals is the default in this study, but some sensitivities to this choice are discussed.Irrespective of the Lagrangian sampling of the points, the final analysis aggregates all points up-and downslope spatially across trajectories.If summarising many cases, the diagnostic  is not averaged over the individual cases but is calculated from the composite of all cases in phase space.In this way, individual sets are weighted with the number of their relevant points.For visualization, points are summarised in histograms but the discretization does not impact the calculation of the topographic forcing diagnostic.

Flow regimes
Two flow indices are used to connect the topographic forcing results to the large-scale circulation.The National Oceanic and Atmospheric Administration (NOAA) provides a daily Pacific-North American (PNA) index based on the National Centers for Environmental Prediction/-National Center for Atmospheric Research (NCEP/NCAR) reanalysis (Kalnay et al., 1996) as part of the Climatic Prediction Center (CPC) (http://www.cpc.ncep.noaa.gov).
The index time series is obtained by applying rotated empirical orthogonal function (EOF) analysis (Barnston and Livezey, 1987) to the monthly mean 500 hPa geopotential height field between 1950 and 2000.The PNA pattern is then identified with the second rotated EOF, which means that the pattern explains a large part of the variability in the northern hemispheric winter atmosphere (Wallace and Gutzler, 1981).In addition, a longitudinally resolved blocking time series based on ERA-Interim is used over Europe.The index is one-dimensional and picks up blocking at the latitude of the climatological storm track.For more details on the blocking index and data, see Hartung et al. (2017).
We composite topographic forcing based on the PNA index, which is representative of the upstream variability, and European blocking is used as an indicator for downstream regimes.Both indices are averaged over the 4 days of each trajectory set and then separated into regimes by their variability (standard deviation).It should be noted that this analysis does not evaluate when during the four-day period the interaction with the topography actually occurs.Regime thresholds need to be surpassed by both the index's average value and its range of variability (standard deviation).

Individual cases
To see how the methodology plays out for specific cases, Although trajectories are slower near the surface, the flow is near-barotropic at higher levels (not shown).The distribution of data in the phase space spanned by the stretching term and the surface forcing (Figure 5a) has an odds ratio indicative of the idealised topographic forcing signal, with an enhanced signal upslope ( u = 8.1,  d = 3.6).It should be stressed again that only points reaching the topography are considered for the calculation of the topographic forcing.Low-level blocked conditions are thus excluded from our analysis, although they are likely to be strongly affected by the presence of topography.This emphasizes that our diagnostic only captures the direct topographic forcing evident in the stretching term in the vorticity equation.January 17, 1982 presents an example of interaction during the negative PNA phase.The large-scale circulation (not shown) is mainly zonal in the free troposphere.Over the 4-day period, the upstream advection changes to be more southwesterly as a trough is extending south over the main mountain ridge.Except at the southernmost latitudes the trajectories converge meridionally and eastward at all levels, crossing the Rocky Mountains around 40-50 • N, and mainly continue downstream at higher levels (Figure 4c,d).The odds ratio of the stretching term is again larger upslope ( u = 3.7) than downslope ( d = 1.4), though in this case the signal of the topographic forcing is reduced (Figure 5b) and  d is close to unity.This indicates that additional processes influence vorticity evolution on the downslope side of the mountain.
The analysis of all 546 trajectory sets includes several examples of atmospheric flow passing straight over the North American Cordillera or meridionally much further north.Besides these cases and variations intermediate between them, several instances of low-level blocked flow occurred (not shown).Thus, air does not cross the mountain ridge in one particular way but rather shows a range of interactions.As the North American Cordillera is latitudinally broad, most trajectories progressing east hit the topography somewhere, although at different elevations.Interactions with topography appear to be, to a large degree, modulated by upstream conditions, which is investigated further in Section 3.3.

Topographic forcing climatology
Results from all trajectory calculations can be summarised in the two-dimensional phase space of the surface forcing and stretching term (Figure 6).Although the shape of the distribution by itself does not clearly suggest a topographic forcing signal, the population numbers are asymmetric These results indicate that the Rocky Mountains induce vortex stretching as expected from theory, but that other processes influence the stretching term as well.Deviations from the assumptions of the idealised theory, like the presence of diabatic and nonlinear processes, are but one example.The impact of choosing Lagrangian sampling, which was motivated by the idealised theory, can be tested by calculating the topographic forcing based on the full three-dimensional fields of divergence and surface forcing.In this case results are in the range of  = 1.3-1.5, which indicates that the method of initial sampling of data matters for the overall result.
Note that the signal of topographic forcing is not as clearly evident if one uses a continuous measure of correlation between stretching term and surface forcing such as the Pearson correlation coefficient.This yields −0.33, which increases to −0.48 if both selection thresholds are increased to 1.5 times the standard deviation.These correspond to coefficients of determination of only 0.11 and 0.23 respectively.
In the left-hand side of Table 1 the results are summarised conditioned on initial latitude, dividing the total range into sets of 5 • lat.Both the values of  and the shape (not shown) of the distribution vary along the latitude composites.The three northernmost intervals exhibit the largest signal of topographic forcing.This increased signal is linked to an elongated shape of the distribution in the direction of the stretching term in phase space (not shown).Farther south, the topographic forcing signal is weaker, although it is still discernible in that the values exceed unity.The ratio of up-and downslope values is also generally decreasing equatorward.An increase of the slope in phase space, and thus of the signal of topographic forcing, with higher geographical latitude is expected from Equation 5 via the Coriolis parameter.Linearly extrapolating an assumed value of  u = 2.8 at 45 • N explains a variability of  u = 2-3.3over the latitude range considered in Table 1.The latitudinal dependence of the topographic forcing signal is thus in line with the underlying simplified theory.Sampling trajectory points based on wind speed (right-hand side of Table 1) provides another perspective.For all cases with wind speed exceeding 10 m⋅s −1 the topographic forcing signal is substantial, and the upslope values consistently exceed the downslope values.For cases with weak winds,  is mostly reduced, although still always exceeds unity.Note that if the threshold on stretching values is eliminated, the values of  are reduced (as expected from the greater amount of noise), but the same dependencies as seen in Table 1 remain (not shown).Further sensitivities to parameters of the method are discussed in Appendix B.

Topographic forcing and flow conditioning
The PNA index captures large-scale circulation variability in the North Pacific.In the case of a positive (negative) PNA the upstream flow is more meridional (zonal) than the climatological wind (Wallace and Gutzler, 1981).However, the North American Cordillera is meridionally broad and in most cases interaction with topography occurs somewhere (see Figure 4).Figure 7 shows the mean topographic forcing and its variability according to a bootstrapping analysis for different phases of the PNA (top row).Considered are 126 positive PNA and 36 negative PNA cases.The number of occurrences in each quadrant q 1 to q 4 (see Figure 2) is bootstrapped 1,999 times while keeping the original sample size, for example, 126 × 1999 positive PNA cases. u, d are then calculated on the re-sampled quadrant population.The upslope topographic forcing signal is overall enhanced during periods of positive PNA ( u = 3.1) compared to negative PNA phases ( u = 2.6) and those of weak PNA pattern ( u = 2.5).With a confidence interval of 95% the difference between the bootstrapped sets of positive and weak PNA occurrences is statistically significantly different from zero (p ≤ .01)upslope of the topography.The difference between positive and negative PNA occurrences is non-zero with p = .06.Downslope, all PNA phase regimes are well separated (p ≤ .01).Again the mean topographic forcing signal is largest during positive PNA phases ( d = 2.1) and reduced during weak ( d = 1.9) and further during negative PNA periods ( d = 1.6).
As the PNA is a statistical measure based on the large-scale pressure pattern, the horizontal wind direction provides a similar criterion for sampling upstream conditions.With a ratio of v/u = 1.12, which compares to 48 • relative to the zonal direction, the flow is orthogonal to the main ridge of the North American Cordillera (fitted based on ERA-Interim topography to be f [] = −49 − 0.89).Figure 8 shows sampling in a reverse fashion to the results in Figure 7, that is, the structure of the flow, using the upstream wind angle, is composited based on .The distribution in dark green (dotted line) represents the initial angle of wind from all trajectory sets at 750 hPa, with a mean angle  = 33 • relative to eastward zonal flow.Restricting the analysis to cases with a strong signal of up-and downslope topographic forcing ( u, d > 2.0, blue dashed line) the mean angle is increased to  = 39 • .The change is due to an overall reduction of negative angles (southward flow) and an increase in northward tilted flow, mainly northeast.Sampling trajectory sets where at least one of the topographic forcing diagnostics is less than one (light green full line), and thus exhibiting a reversal of the expected distribution, indicates increased occurrence of cases with a southward orientation of the initial flow and a mean angle of  = 13 • .Dividing the composites based on only  u or  d alone points to the origin of the overall division mainly coming from downslope regions (not shown).This is in line with the PNA composites producing larger variability of  among regimes downslope than upslope.
Figure 7 (bottom row) also shows the topographic forcing signal sampled for European blocking where the interaction with the topography is leading the downstream regime by 1 day.The upslope topographic signal forcing is effectively indistinguishable between periods of enhanced and reduced blocking ( u = 2.6, p = .39)but  d is increased during periods of enhanced blocking from  d = 1.9 to  d = 2.1 (p = .02).In the range of lags between 0 and 4 days (not shown), the downslope difference between periods of strong and weak blocking is largest for a lag of 1 day and otherwise the p-value is exceeding the 0.05 threshold.In contrast, the difference in the signal of upslope topographic forcing peaks 2 days before periods of enhanced blocking (p = .06)and is only slightly reduced for lags of 3 or 4 days (p = .24).

DISCUSSION, CONCLUSIONS AND OUTLOOK
The main findings from this study are that the large-scale topographic forcing signal (quantified as the odds ratio of the expected sign of vortex stretching from the topographic forcing) associated with the North American Cordillera is (a) evident in reanalysis data, (b) typically larger upslope than downslope, and is increased when trajectories are (c) located further north or (d) associated with larger wind speeds.In addition, (e) upstream conditions of positive PNA increase the signal of topographic forcing and (f) the largest signal of topographic forcing is associated with northeast tilted flow.Finally, (g) more frequent European blocking is associated with stronger topographic forcing downslope of the topography.A variety of overflow patterns exist and from this analysis it cannot be concluded that air flows more often over or around the North American Cordillera.
The goal of the analysis was to determine how large the signal of topographic forcing is in an atmospheric dataset.Although a signal is present and its characteristics are physically consistent, the signal is not as clear-cut as in the idealised model.Possible reasons are deviations from the idealised theory (non-zonal upstream and non-single layerflow, barotropic potential vorticity not conserved e.g.due to baroclinic processes, … ) and the fact that the surface forcing is a local parameter and thus does not include potential spatial delays of the topographic forcing signal at higher levels.Given all these potential complications, it is perhaps surprising (but also reassuring) that the textbook theory appears to describe the gross aspects of the flow behaviour.Although the ERA-Interim reanalysis provides an approximation to reality, it is likely still subject to errors of the numerical formulations in the underlying IFS atmospheric model.In particular, limitations in the lower-atmospheric drag parametrizations (Sandu et al., 2016) can influence the results because low-level winds are not assimilated over the continents (Dee et al., 2011).The origin of the different characteristics of  u and  d has not been investigated in detail.Considering the formulation of the linearised boundary condition w ′ = u h ′ x (Valdes and Hoskins, 1991) and the shape of the North American Cordillera it can be noted that the gradients are typically larger upslope than downslope (not shown).This feature is partly due to the sharp coastal boundary to the west in contrast to the slope extending into the main continent to the east.It is furthermore expected that reduced surface drag over the upstream ocean, compared to over land, causes increased upslope winds.These observations together indicate that the enhanced signal of upslope topographic forcing is likely caused by stronger surface forcing.Another explanation could be the development of flow parallel to the topography downslope of the topography (Hsu, 1987), which shields the westerly flow at higher levels from the surface forcing and prevents downslope topographic stretching.The lack of descent back to the original level is for example evident in the lee of the North American Cordillera in Figure 4c.Notwithstanding the difference in magnitude between  u and  d , they both exhibit a sensitivity to the initial latitude as predicted by the simplified model of the topographic forcing (Equation 5), which provides a physical plausibility to our diagnostic.
Sensitivity of the topographic forcing signal to upstream conditions, such as PNA phase and wind angle, provides a mechanism for remote impacts of upstream stationary wave patterns (e.g.Hoskins and Karoly, 1981).The atmospheric conditions upstream of the North American Cordillera have been shown to be influenced by diabatic heating of the Pacific storm track (Valdes and Hoskins, 1989), and Saulière et al. (2012) show that adding the Tibetan Plateau reduces the downstream response of the Rocky Mountains in an idealised model.Downstream results point to a link between an increased downslope signal of topographic forcing and an enhanced occurrence of atmospheric blocking shortly thereafter.With a Rossby wave group velocity of 2u = 50 m⋅s −1 the influence of the topography travels about 4,500 km (60 • lon) downstream per day, and thus the dynamic link between the North American Cordillera and European climate on the identified time-scale is physically plausible.It should also be noted that in the current set-up, where blocking indices are aggregated over four-day periods, it cannot be determined on which day the interaction with the topography occurs.The true lag thus might be different from the one presented here.However, studies by Pithan et al. (2016) and van Niekerk et al. (2017) also point to a link between the representation of topography and downstream flow conditions.Pithan et al. (2016) detect more European blocking when the low-level drag in the Unified Model is increased.Increased low-level drag presumably increases the topographic forcing signal, although these results are not specifically linked to the North American Cordillera.van Niekerk et al. (2017) similarly find a more tilted North Atlantic jet with stronger low-level blocking in the Canadian Atmospheric General Circulation Model (AGCM) CanAM.
The method presented here to diagnose a topographic forcing signal can also be applied in a numerical weather prediction or climate model.As a climate model likely produces a different climatology of circulation regimes it would be necessary to sub-sample the topographic forcing signal based on up-and downstream conditions to ensure a meaningful comparison.The topographic forcing diagnostic can be used to highlight model biases and sensitivity to model resolution and parametrizations.Considering a general comparison among models it might be necessary to further develop the method to produce similar results from three-dimensional Eulerian sampling as with sampling along trajectories.Such an extended method would need to ensure that air is actually crossing and interacting with the topography and exclude stagnant parts of the flow from the analysis.The study presented here suggests that such an effort would be beneficial for studying model biases related to the representation of topography.

APPENDIX A: Shallow-water model
We use a two-dimensional shallow-water model (Hogan, 2018) with three prognostic variables: the depth of the fluid layer H and the horizontal wind velocities u and v.In a channel with cyclic x-boundary conditions, mean westerlies (u = 20 m⋅s −1 ) flow over a localised Gaussian ridge h with maximum height h max = 2000 m: The model grid has 260 zonal and 60 meridional points with dx = dy = 10 5 m.The meridional extent of the idealised topography is chosen to roughly match the extent of the North American Cordillera.The width is relatively increased to allow more samples over elevated areas.Conclusions drawn from the results do not depend on the shape of the idealised ridge assumed in Equation A1 (not shown).The central latitude is 45 • N and the depth of the flow is defined as H = H 0 + h top − h with the initial depth H(t = 0) = H 0 = 10 4 m and the depth of the flow exceeding H 0 -h denoted as h top .The shallow-water equations on the -plane are solved with a Lax-Wendroff scheme and a time step of Δt = 3 s over a simulation time of 7 days.This means that the Rossby wave initiated by topography can move about half-way around the channel so that the wind directly upstream of the topography remains zonal.At the analysed time the flow is in quasi-steady-state above the topography.

APPENDIX B: Parameter sensitivities
The proposed approach to diagnose the signal of topographic forcing via the stretching term is potentially sensitive to subjective parameter choices.Trajectories can be restricted vertically, both reducing the maximum height (default 300 hPa), and limiting the number of layers from which trajectories are initialised (default between 850 and 600 hPa).Increasing the threshold of topographic height over which points are considered (default is all points, threshold 0 m) limits the analysis spatially to the main mountain ridge.
If all trajectory points are taken into account, even those associated with weak stretching and small surface forcing, the impact of most of the above-mentioned Figures 4 and 5 show examples of individual trajectory sets for two phases of the PNA.From 1 January 1981 the evolution of the large-scale circulation (not shown) is dominated by a ridge over the central North American Cordillera during a period of positive PNA.The pressure ridge is associated with southerly/southwesterly flow upstream of the main mountain ridge which is driving the circulation for most of the period.This flow pattern only breaks up during the course of the last day of the 4-day trajectory period.Trajectories indicate meridional flow upstream of the coastline and most of the interaction with the terrain occurs north of 60 • N (Figure 4a,b).

F
Atmospheric trajectories initialised at (a,c) 850 hPa and (b,d) 650 hPa for release dates (a,b) 1 January 1981 and (c,d) 17 January 1982.Colours denote the pressure at each point.Topography in grey contours every 500 m F I G U R E 5 Phase space of stretching term versus surface forcing for (a) 1 January 1981 and (b) 17 January 1982.Colours indicate the density of points.Numbers in the corners denote the total number of observed points in each quadrant outside of the grey bars, which represent the thresholds below which data are excluded from the analysis

F
I G U R E 6 Summary of 546 winter cases in phase space of stretching and surface forcing term.Colours indicate the density of points.Numbers in the corners denote the total number of observed points in each quadrant outside of the grey bars, which represent the thresholds below which data is excluded from the analysis both up-and downslope.Disregarding the smallest surface forcing values gives odds ratios of  u = 2.7 and  d = 2.0.Thus, upslope compression is more clearly detected and downslope values are reduced by about 25% in comparison.

F
Results of the odds ratio  from 1,999 resamples of data, conditioning with respect to flow regimes PNA and European atmospheric blocking.Boxes show interquartile range, whiskers show the range of 95% of all data.PNA phases are distinguished as positive (+), negative (−) or weak regime pattern (0).Thresholds between regimes are determined based on one standard deviation of variability F I G U R E 8 Normalised distribution of initial wind angle at 750 hPa.The distribution of all 546 cases is shown by a dotted dark green line.Trajectory sets with both  u and  d larger than 2 are shown in a dashed blue line, and those with at least one diagnostic smaller than 1 in a full light green line.The short vertical lines at the bottom of the figure indicate the mean value of the respective distribution

Lat [ • N] 𝝉 u 𝝉 d |v| [m⋅s −1 ] 𝝉 u 𝝉 d
Topographic forcing signal from combination of all cases for subsets depending on initial latitude and total wind speed.The mean values are:  u = 2.7,  d = 2.0 TA B L E 1