C-Band Dual-Doppler Retrievals in Complex Terrain: Improving the Knowledge of Severe Storm Dynamics in Catalonia

: Convective activity in Catalonia (northeastern Spain) mainly occurs during summer and autumn, with severe weather occurring 33 days per year on average. In some cases, the storms have unexpected propagation characteristics, likely due to a combination of the complex topography and the thunderstorms’ propagation mechanisms. Partly due to the local nature of the events, numerical weather prediction models are not able to accurately nowcast the complex mesoscale mechanisms (i


Introduction
Severe weather (defined for the area of study as large hail with diameter over 2 cm, straight-line winds greater than 25 m•s −1 associated with convective thunderstorms, downbursts, tornadoes, and/or waterspouts, [1]) is quite frequent in Catalonia, a small region of 32,000 km 2 in the northeast of the Iberian Peninsula (see, for instance, [2,3]).This area has complex topography: the Pyrenees to the north reaching an altitude higher than 3000 m, the Central Depression to the west, and the Littoral and Pre-Littoral ranges to the east, which all play important roles in the development and evolution of these severe thunderstorms [4,5].Although the topography and proximity to the sea favors the formation of such interesting weather, it also represents a handicap in surveillance tasks, since operational numerical weather prediction models have many difficulties in reproducing the topography.This does not only affect the location and intensity of thunderstorms in convection-allowing models but also the capability to reproduce their fast-changing thunderstorm-scale-related interactions, which represent a major factor in nowcasting severe weather.Because of this, remote sensing, especially weather radar, serves as a fundamental tool for observing complex severe weather events in areas with complex terrain and highly variable atmospheric conditions [6].The complexity of the Catalan orography and the limitations in tools establish that the forecast temporal limit is two hours, for which a severe weather phenomenon is predicted and located with acceptable accuracy.This means that thunderstorm-track probability cones from radar echo extrapolation are retrieved out to 2 h, although the accuracy degrades after 30 min.
An important technique for studying thunderstorms, which has also become a pursuit in several weather services worldwide, is the multiple-Doppler radar retrieval-obtaining 3D winds within a thunderstorm from radar observations.Retrievals of the entire wind field within a storm allow an understanding of how it may interact not just with its environment but also with its own internal evolution during its life cycle, and therefore, allow for the potential to improve short-term forecasting skill.Furthermore, the derived fields may be used as input for high-resolution models, which can allow us to forecast and understand complex cases in the future [7,8].Deriving the 3D wind field with multiple-Doppler techniques can be achieved using the observed radial velocity from two or more non-collinearly located Doppler radars.Different techniques using cylindrical and Cartesian coordinates system are explained, for instance, in [9,10], which are the basis for more advanced techniques presented in the following paragraphs.The mean radial velocity from every radar is related to the rectangular components of the mean velocity of the hydrometeors (u, v, w) through the azimuth and elevation angles.Therefore, it is possible to obtain a unique solution for the horizontal and vertical velocities at a point by combining the radial velocities in a closed equation system.
Usually, due to the availability of radars in research projects, one of the most used multiple-Doppler techniques is based on the observations of only two Doppler radars; this is the dual-Doppler technique.Since this is also the focus of the present work, hereafter the paper is focused on this dual-Doppler technique.In the dual-Doppler configuration, and because only two weather radars are used to retrieve three components of the wind field, a major handicap arises when computing the solution from the system of equations; it becomes under-determined (two main equations and three unknowns).This has to be solved by applying boundary conditions on the vertical velocity and by imposing the mass continuity equation, making assumptions on the relationship between the vertical velocity and the terminal velocity of the particles [11].Furthermore, other features hamper the easy obtention of the wind field.It is well known that weather radar products may be affected by sampling errors, beam effects, blocking, changes in spatial resolution when interpolating the data, among others [12,13].These effects may impact the retrieval of the integrated divergence and the vertical component of the velocity, where most of the assumptions and boundary conditions are applied.Therefore, there is a need to find a technique that eases such a process.One of the techniques that has been shown to lessen the handicaps mentioned is the variational method, which consists of determining the solution by minimizing a cost function similar to a data assimilation process (i.e., [14][15][16]).In the variational method, wind retrieval is solved globally by minimizing the difference between the observed radial velocities and the analyzed wind solution, with additional constraints from mass continuity and spatial filtering.Since the variational method has been demonstrated to give better results than the traditional solution [16][17][18] and the computational resources are increasing in research centers, it has become one of the most used methods today (see, for instance, [19,20]).
As discussed, dual-Doppler radar techniques for obtaining wind fields have been continuously developed and enhanced, facilitating the obtention of important results such as the dynamics of supercells, tornadoes, or hurricanes [21][22][23][24].Those results have been mostly obtained during field campaigns by means of transportable and mobile research radars.Therefore, there are just a few studies testing the techniques with operational systems, which face some difficulties for retrieving 3D wind fields.For instance, dealing with the "Doppler Dilemma" in an operational arena might force scanning in a longer velocity range to avoid aliasing of radial velocities outside the Nyquist interval [25,26], thus, having a shorter unambiguous range for surveillance tasks, and vice versa.Additionally, radar baselines (separation between radars) are usually around 100 km or more, which impacts the coverage area and the data resolution [27].Furthermore, in some regions, including the area of study, operational systems have to deal with complex topography that can compromise the data collected due to beam blockage, or even the possibility (rare but plausible) that the scanning strategies may not be synced in time.These are just some of the difficulties that dual-Doppler techniques encounter when retrieving the 3D wind field with operational radars.
A near-real-time dual-Doppler analysis was proposed for the first time in [28], implementing the technique in the Mesoscale Alpine Program (MAP) field campaign in 1999.This technique was based on a complex terrain-adapted version of the Multiple-Doppler Synthesis and Continuity Adjustment Technique (MUSCAT, [29]), with a radar baseline of about 70 km and two C-band radars.Later, the MUSCAT algorithm was demonstrated successfully to produce 3D wind syntheses and it was temporarily implemented for the French operational radar network [30].A near-real-time integrated display and analysis tool was also proposed in [31] in order to help scientists in field projects and nowcasting in warm-season convection, using dual-Doppler techniques in Colorado, USA, with four radars, from which two of them were operational.Furthermore, the capability to retrieve multiple-Doppler winds was also demonstrated by using the monostatic operational German Weather Service's (DWD) radar network and a bistatic radar network [32] as a possible operational surveillance tool around the Frankfurt/Main International Airport [33].
The advances in recent years in the use and development of near-real-time multiple-Doppler techniques with operational radar networks, as well as the improvements achieved in complex terrain, have provided the motivation and foundation for the study presented here.The main purpose of this paper is to demonstrate how the operational radar network, XRAD (Xarxa de RADars, in Catalan) of the Meteorological Service of Catalonia, hereinafter SMC, is suitable to retrieve dual-Doppler winds by using a free variational Doppler software within the Lidar Radar Open Software Environment (LROSE) project [34]: the Spline Analysis at Mesoscale Utilizing Radar and Aircraft Instrumentation (SAMURAI).We will demonstrate the capability of this tool to retrieve winds in a complex terrain area and how it represents a big step at the SMC, being the first service to retrieve dual-Doppler analysis at the Iberian Peninsula and to obtain dynamical features within particular severe storms with unexpected motion patterns that have affected the Catalan area.This manuscript is divided as follows: the motivation of the present work is explained in Section 2. Section 3 describes the data and methods used: the XRAD network and data and complexity of the area of study and the SAMURAI technique.Section 4 shows the results for two severe weather events affecting the area-a splitting severe hailstorm in July 2013, and a tornado-producing supercell channeling through a valley in January 2018.Finally, the discussion and conclusions are given in Section 5, where future work is also proposed.

Motivation
Severe weather phenomena in Catalonia, especially hail-related events, cause significant economic losses in infrastructure and agricultural yields, and they occur with high frequency.For instance, the authors in [35] found that the yearly average economical losses produced only by hail events in the western area of Catalonia were EUR 15M for the period 2000-2009.In the case of straight-line winds or tornadoes, the frequency of occurrence is lower, although most of the events affect densely populated areas and therefore, the economic impact may be even higher [36,37].
Figure 1a shows the number of severe weather events per km 2 recorded per county for the period 2001-2019.This map was constructed using an internal database of the SMC.The SMC database includes a total of 275 events that have been recorded within the period, where around 67% of the cases are hailstorms and the remaining 33% are wind-related phenomena (14% downburst and wind gusts and 19% tornadoes).The severe storm climatology is primarily obtained from storm spotters' observations and this requires consideration of several factors when analyzing the database.Such considerations include the dependency of the spotters' observations on population density and the reliance and vulnerability of social networks [38,39].In particular, the latter has favored the increase in the number of records during the most recent years: the increase in social media platforms and the ease to share live videos and pictures has resulted in 50% of the total events being recorded within the last six years.
western region represents a hotspot.This is mainly because it coincides with a major agricultural area that has been the target for several research and stakeholders' projects [1,35], providing the area with a hail pad network and with citizen science campaigns to enhance the warning systems.Finally, it is worth mentioning another area that is clearly underrepresented within the severe weather database: the central and northern area, which is the mountainous region of Catalonia (Pre-Pyrenees and Pyrenees).From the previous analysis, it might seem that the region has a low frequency of occurrence of severe phenomena, although a large number of the most recent events have occurred in this area.Besides, some cases which are currently under analysis seem to show that the complex topography of the region is not preventing the creation or propagation of severe storms but, on the contrary, it enhances their longevity or even helps to present anomalous patterns of propagation.This is the case, for instance, of two events that are presented in the Results section: a splitting storm on 10 July 2013 and a tornado-producing supercell that was channeled through a valley on 7 January 2018.In respect to the nowcasting techniques at the SMC, there are two efforts working in parallel.The first one is based on mesoscale numerical weather prediction (NWP) models: the Weather Research and Forecasting model (WRF, [40]), run at 3 km resolution, which is the finest operational resolution at the SMC.During this process, data from METeorological Aerodrome Reports (METARs, from airports), automatic weather stations (AWS), and the weather radar are assimilated into the simulations, resulting in improved forecast skill [41,42] for events related to meso-alpha systems and We can distinguish three main hotspot severe weather areas in Catalonia, represented in Figure 1a.The regions with the highest number of phenomena are the internal western region (the Plains), the central coast, and the northeastern counties.These areas have counties with more than 0.015 events•km −2 (in maroon in Figure 1a).This is mainly due to hail reports in the western plains and in the north-central counties.On the contrary, wind-related phenomena have been observed mainly along the coastline and in the northeastern counties.The severe weather reports' spatial distribution may be explained by a higher number of occurrences and a higher density of population (Figure 1b) reporting impacts in infrastructure and agricultural yields.That is, people report more cases when there are personal economic losses.Having said that, we should expect a lower occurrence in the western and northwestern areas, since this also coincides with a sparsely populated area.The central western region represents a hotspot.This is mainly because it coincides with a major agricultural area that has been the target for several research and stakeholders' projects [1,35], providing the area with a hail pad network and with citizen science campaigns to enhance the warning systems.Finally, it is worth mentioning another area that is clearly underrepresented within the severe weather database: the central and northern area, which is the mountainous region of Catalonia (Pre-Pyrenees and Pyrenees).From the previous analysis, it might seem that the region has a low frequency of occurrence of severe phenomena, although a large number of the most recent events have occurred in this area.Besides, some cases which are currently under analysis seem to show that the complex topography of the region is not preventing the creation or propagation of severe storms but, on the contrary, it enhances their longevity or even helps to present anomalous patterns of propagation.This is the case, for instance, of two events that are presented in the Results section: a splitting storm on 10 July 2013 and a tornado-producing supercell that was channeled through a valley on 7 January 2018.
In respect to the nowcasting techniques at the SMC, there are two efforts working in parallel.The first one is based on mesoscale numerical weather prediction (NWP) models: the Weather Research and Forecasting model (WRF, [40]), run at 3 km resolution, which is the finest operational resolution at the SMC.During this process, data from METeorological Aerodrome Reports (METARs, from airports), automatic weather stations (AWS), and the weather radar are assimilated into the simulations, resulting in improved forecast skill [41,42] for events related to meso-alpha systems and therefore, provide a valuable tool for surveillance purposes.The second procedure considers exclusively remote sensing data sources and runs operationally every six minutes.This combines lightning jump (LJ) warnings [1] with an object-based identification, characterization, tracking, and nowcasting algorithm of thunderstorms to analyze the volumetric data of the radar network [6,43] to generate a cone of probabilities of occurrence of severe weather in the next 2 h.

The XRAD Network and Data
The radar network of the SMC, the XRAD, is composed of four C-Band (wavelength of 5.36 cm) Doppler radar with single polarization (see [44], for more details).Table 1 shows the main features of the XRAD, while Figure 2a presents the spatial distribution and coverage.The complex topography of the territory causes some partial or total blockage of the lower elevation beams for individual radars [45] (areas in black in Figure 2).Since the main purpose of the XRAD is to provide the information necessary for surveillance tasks during adverse weather episodes, the network was configured to provide a good reflectivity composite with almost no blockage under 3000 m ASL (approximately the top of the Pyrenees).In spite of the processor filters, there remain some further corrections that must be made before ingesting the data into the dual-Doppler algorithm.In order to overcome the aliasing issue, that is, when the scatterers move faster than the Nyquist velocity and the estimated velocity is folded back to fit in the Nyquist velocity range [25], the commercial radar processor applies the dual-PRF technique [47].This technique also causes problems in cases of high azimuthal shear, quite common in severe weather events, resulting in velocity outliers and errors [48].For this reason, the radial velocity data are processed with a new algorithm to correct such errors when the dual-PRF technique The operational chain, managed using Sigmet-Vaïsala IRIS© software [46], retrieves each single radar volume in a 6-min cycle, composed of three different tasks: VOL_A, VOL_B, and VOL_C.
Since the VOL_A is a long-range surveillance task which stores only reflectivity, the radar data used in the work presented here are directly from VOL_B and VOL_C volumes, which also include the radial velocity.The reflectivity and radial velocity fields are automatically filtered in the processor during a real-time quality control procedure.This allows us to obtain the best and fastest field possible to combine and develop the surveillance composite, removing the main non-meteorological targets.The automated filter applies Signal-to-Noise Ratio (SNR), Clutter-to-Signal Ratio (CSR), and Signal Quality Index (SQI) [46] thresholds to each of the range bins in the radar data after spectral moment estimation, discarding those bins below the thresholds.These thresholds were set manually by the radar technicians at the SMC and vary slightly depending on the radar and the task (see [44] for a detailed explanation of the values used in the thresholding).However, and especially in high azimuthal shear situations, these threshold values result in a very strict data removal, especially for the radial velocity, which then causes missing bins in the radar data (an example is provided in the Results section).Furthermore, neither the raw unfiltered nor the spectral width (distribution of velocities in a radar pixel) are archived.These missing data present a handicap to overcome with the dual-Doppler analysis.
In spite of the processor filters, there remain some further corrections that must be made before ingesting the data into the dual-Doppler algorithm.In order to overcome the aliasing issue, that is, when the scatterers move faster than the Nyquist velocity and the estimated velocity is folded back to fit in the Nyquist velocity range [25], the commercial radar processor applies the dual-PRF technique [47].This technique also causes problems in cases of high azimuthal shear, quite common in severe weather events, resulting in velocity outliers and errors [48].For this reason, the radial velocity data are processed with a new algorithm to correct such errors when the dual-PRF technique is applied [26,44].Finally, before starting the dual-Doppler analysis, data are manually edited using SOLO software [49].This is done to manually correct errors such as ground clutter, second trip echoes, or electromagnetic interferences [25], especially around the area of the target storm being analyzed.

Dual-Doppler Configuration
According to [13], in order to resolve storm-relative winds and updrafts, the paired radars must be able to resolve horizontal scales finer than 3-6 km.This is affected by the distance between the paired radars (baseline): if the radars are too far apart, the data will be degraded at the edges of the dual lobes because of beam spread and elevation.Furthermore, to obtain the optimal dual-Doppler configuration and due to the maximum acceptable error in the velocity (setting an upper limit of 3 m•s −1 in the standard deviation of wind velocity [33,50]) as well as the total area that we want to cover, the best beam-crossing angle between the radars must be chosen [27].In this sense, taking into account that the XRAD radars' beam width is 1.2 • and that we cannot change the operational scanning strategies, in the present study, the beam-crossing angle is set to 30 • and is used for the following paired radars: CDV-PBE and CDV-LMI (Figure 2b).
Table 2 summarizes the dual-Doppler coverage, baseline, gate resolution of the farthest point (less than or equal to 3 km for both pairs), and blockage for each pair of radars (chosen ones are shaded in grey).The total area (Table 2) indicates the maximum coverage within the dual lobes.However, since the convective storms in Catalonia are very localized phenomena, we aimed to set the desirable gate resolution (the maximum permitted) to 1.7 km, which allows us to cover optimal areas of 10,340 and 7994 km 2 for the CDV-PBE and CDV-LMI pairs, respectively.This coverage allows us to retrieve data with good confidence in the regions more affected by severe weather in Catalonia as we saw in Figure 1: the western Plains, the central area, and the central and southern coast (Figure 2b).However, this leaves the northeastern area uncovered due to the large baseline from the PDA radar to the closest radars (109.01 km to PBE and 136.10 km to PDA, as seen in Table 2).The pair PBE-LMI was not chosen due to the overlap of coverage with the chosen pairs, with a configuration resulting in a lower resolution of the data.Finally, the topography blockage for the lowest beam is also low in both chosen pairs (Table 2), since most of the blocked area falls within the blind area between the dual lobes (Figure 2b), where the wind cannot be resolved correctly due to the collinearity of the radial wind components of each radar.This is less than 2% of the total area in both cases, and less than 1.5% in the CDV-PBE pair.

SAMURAI
The SAMURAI software [51] is a three-dimensional variational analysis method which is able to retrieve kinematic and thermodynamic features by ingesting different meteorological data sources.The algorithm determines the most probable state of the atmosphere based on observations by minimizing a variational cost function in incremental form and using cubic B-spline finite elements [52][53][54].This method can incorporate complex and multiple sources of observations that allow it to have a better analysis and background estimates (e.g., ground-based and aircraft radar data, atmospheric vertical profiles from soundings, NWP model fields, and AWS data).In SAMURAI, the analyzed wind field is represented as a function composed of the cubic B-spline finite elements rather than a grid point representation [53].Since the basic functions are cubic polynomials, the spatial gradients of the wind field can be obtained analytically.The direct specification of the gradients via spline representation of the wind field allows the mass continuity constraint to be enforced without explicit integration through the use of "pseudo-observations" in the domain interior.Through this technique, mass continuity is satisfied everywhere in the domain without requiring integration from the lower or upper boundary condition.Since the solution is solved globally through the variational constraint, the magnitude of the retrieved vertical velocity is still influenced by divergence throughout each column, but the spline-based approach minimizes the impact of the boundary conditions.Further details on SAMURAI can be found in [54,55].
SAMURAI utilizes low-pass filters to reduce noise in the solution and also provide some limited interpolation in areas of missing data to produce a smooth, continuous wind field.The amount of filtering is specified by the user through the prescribed filter lengths of a spline 'cut-off' that minimizes the third derivative of the wind field [56] and a Gaussian low-pass filter [57].It is important to set the filter parameters correctly, since high values may result in an over-smoothed analysis and loss of fine scale features, while insufficient filtering can result in a noisy analysis [27,54,55].The prescribed values for the analysis of the study cases are presented in Table 3.
Table 3. Spline Analysis at Mesoscale Utilizing Radar and Aircraft Instrumentation (SAMURAI) main configuration for the study.∆ is the length scale used for the filters (horizontal and vertical grid scale).

Z Radius Infl. (km)
Gaussian Filter (∆) Spline Cutoff (∆) Analysis Grid Increment (km) 1 4; 4; 2 2; 2; 2 0.5; 0.5; 0.5 Since the current version of SAMURAI does not enforce the 30 • beam-crossing angle for the dual-Doppler automatically, this was enforced by picking the analysis cases that match automatically this requirement, so the storm or system falls within the required lobes.After a sensitivity test, the Gaussian recursive filter length scale in grid points was set to 4 ∆ in the horizontal directions (X and Y) and 2 ∆ for the vertical (Z), while the spline cutoff length scale, also in grid points, was left to 2 ∆ in all directions, as recommended by [51].The analysis grid spacing was set to 0.5 km, in order to resolve updrafts smaller than 1 km wide.The radius of influence for gridding the reflectivity field was set to 1 and no reflectivity values were masked due to the existence of gaps in the raw data.It is important to note that the resolution shown in Table 2 is the minimum for each radar-this is at the furthest point.In both analyzed cases, the storm was 15-30 km away from the closest radar (CDV), which gave a gate resolution of 0.4-0.6 km.The analysis grid spacing of 0.5 km is consistent with the original gate resolution, and with the prescribed filtering, the analysis resolves features with spatial scales larger than 2 km in the horizontal and 1 km in the vertical.
Since our analysis has been purely based on observed radar data, to test the capabilities of the configuration of SAMURAI and XRAD, the only data used have been the reflectivity and radial wind from the radars selected (no NWP or other data sources were used).
In order to provide a reasonable density field for the mass continuity constraint, a vertical sounding is required for the analysis.The only sounding station in Catalonia is located in Barcelona, very close to the coast, which may have moist maritime influence in low levels.For this reason, for the case located in the western Plains (inland), we have used a WRF forecast sounding of the closest location (shown in the Results section).Sensitivity tests demonstrate that the sounding density profile does not have a significant impact on the analysis (not shown).Since there was no reliable first guess wind field available, the first guess background field wind was set to zero and the background error was set to a high value, so that the wind solution was determined solely by the observations, mass continuity constraint, and filtering.Finally, the storm motion resulting from the tracking system [43] was introduced in the advection scheme, to take into account the mean storm propagation (steering wind) in the analysis process.In addition, data were converted to CfRadial format [58], with the same LROSE package tools, before ingesting it into the SAMURAI software.

Results
The results section is divided into two subsections to illustrate the advantages and capabilities of dual-Doppler retrievals in the Catalan area.First, a tornado-producing supercell case on 7 January 2018 is used to demonstrate the performance of the SAMURAI technique when retrieving storm dynamics over complex terrain.Second, a splitting storm on 10 July 2013 is analyzed to show the possible role of topography associated with an anomalously propagating system.

SAMURAI Performance
An unusual weather system occurred over Catalonia on 6-7 January 2018.The severe weather event occurred outside the convective season, which, in the northwestern Mediterranean, is between May and September, and it occurred overnight, while the usual thunderstorm time of occurrence in Catalonia is between 12 and 19 UTC [3,5,59].Several rain bands with predominant southeasterly flow affected the area.Within those bands, embedded convective cells produced severe weather, causing two significant tornadoes, both rated EF2 (Enhanced Fujita scale [60]).
The state of the atmosphere for the 2018 tornado event is shown with the observed Barcelona sounding at 00 UTC (Figure 3a), 43 km southeast of the first supercell storm location (Figure 3b).This sounding profile is the same used in the SAMURAI analysis for the event.During the 6-7 January 2018 event, a deep low visible at different pressure levels was located in the center of the Iberian Peninsula, placing the exit region jet at 300 hPa over Catalonia, with strong southeasterly flow, as can be seen in Figure 3.The axis of the low was tilted towards northeast in lower levels, positioning the center in the eastern border of Catalonia, which favored a more easterly flow at low levels, causing warm moist advection from the Mediterranean Sea in the boundary layer (Figure 3).The high instability of the atmosphere was demonstrated with a Most Unstable CAPE (MUCAPE) of 1596 J/kg, which was actually the same value found for a surface-based CAPE (SBCAPE).This partially explains why the main rain bands that contained the two EF2 tornadoes originated over the sea and once they touched inland, they organized into supercells.Furthermore, the clockwise curved hodograph at low levels depicted in Figure 3, and the deep layer bulk shear parameters values reaching more than 40 kt demonstrate the warm advection mentioned before and the strong shear necessary for supercell formation [61].This is also shown with the high values of Storm Relative Helicity (SRH) in Figure 3, with more than 100 m 2 s −2 at lower levels, which is sufficient for tornado formation.
Remote Sens. 2019, 11, x FOR PEER REVIEW 9 of 26 which was actually the same value found for a surface-based CAPE (SBCAPE).This partially explains why the main rain bands that contained the two EF2 tornadoes originated over the sea and once they touched inland, they organized into supercells.Furthermore, the clockwise curved hodograph at low levels depicted in Figure 3, and the deep layer bulk shear parameters values reaching more than 40 kt demonstrate the warm advection mentioned before and the strong shear necessary for supercell formation [61].This is also shown with the high values of Storm Relative Helicity (SRH) in Figure 3, with more than 100 m 2 s −2 at lower levels, which is sufficient for tornado formation.One of the tornadoes (red stars in Figure 4) affected the central region of Catalonia, damaging farms and the surrounding pine woods, and the supercell propagation was interesting because it was channeled completely through a narrow valley.This can be seen in Figure 4, which shows the track of the convective cell that produced the tornado, overlaid on top of the topography.Due to the resolution and the scanning strategies, it has never been possible to observe the tornado itself with the XRAD (smaller than 1 km wide), and thus, as in analyzed cases from the past, the tornado information was obtained with "on the ground" post-event debris analysis [37].Additionally, in this case, the storm propagated up the valley (valley floor height of 250 m ASL), and the closest radar, CDV, was located 42 km away at 825 m ASL.This implies that the lowest radar beam was at 440 m AGL with respect to the radar, which is ≈ 1000 m ASL relative to the valley ground level.Although the radar was not able to resolve the tornado vortex signature (TVS), the mesocyclone was visible in several consecutive volumes of the radar data (green dots in Figure 4).One of the tornadoes (red stars in Figure 4) affected the central region of Catalonia, damaging farms and the surrounding pine woods, and the supercell propagation was interesting because it was channeled completely through a narrow valley.This can be seen in Figure 4, which shows the track of the convective cell that produced the tornado, overlaid on top of the topography.Due to the resolution and the scanning strategies, it has never been possible to observe the tornado itself with the XRAD (smaller than 1 km wide), and thus, as in analyzed cases from the past, the tornado information was obtained with "on the ground" post-event debris analysis [37].Additionally, in this case, the storm propagated up the valley (valley floor height of 250 m ASL), and the closest radar, CDV, was located 42 km away at 825 m ASL.This implies that the lowest radar beam was at 440 m AGL with respect to the radar, which is ≈ 1000 m ASL relative to the valley ground level.Although the radar was not able to resolve the tornado vortex signature (TVS), the mesocyclone was visible in several consecutive volumes of the radar data (green dots in Figure 4).Figure 5 shows the raw reflectivity field, where a clearly defined hook echo can be distinguished in different radar volume scans (upper panels).The mesocyclone of the supercell is demonstrated with the radial velocity data (lower panels in Figure 5), where the typical couplet of velocity is depicted.In this case, the raw velocity data are overly thresholded; the filtering of the raw velocity data during the operational processing results in a loss of practically half of the mesocyclone in all the volumes (white pixels in dashed circles).As explained in Section 3.3, optimal tuning of the filters within SAMURAI allows us to fill data gaps during the analysis and therefore, to minimize the effect of the previous stringent operational filtering problem.The retrieved reflectivity and horizontal wind field at 1 km height is depicted in Figure 6a, where the southeasterly steering flow (Figure 3) is shown by the thick black arrow.A large area of weaker reflectivity values surrounds the main supercell, indicating that the supercell was embedded in the precipitation band, and the tornado was mostly rain-wrapped.In addition, Figure 6a depicts a region at the southeast corner where the wind is slightly more intense and converging within the rear part of the supercell, which is consistent with the storm inflow feeding the main updraft region (red shaded area in Figure 6c, indicating positive vertical velocity).The SAMURAI retrievals also resolve the closed circulation in the storm-relative horizontal winds (black arrows in Figure 6b), which depicts the mesocyclone at mid-levels (2.5 km AGL), and the area with a couplet of positive and negative vertical velocity (red and blue in Figure 6c).This vertical velocity couplet depicts the main updraft (red shaded) and the corresponding forward-flank downdraft (FFD; blue shaded) of this rain-wrapped heavy precipitation (HP) supercell [64][65][66].However, Figure 6c also depicts a complicated vertical vorticity retrieval, in 10 −5 s -1 , where the main couplet is displaced from the updraft-downdraft area, and therefore, we cannot completely relate the mesocyclone with the main updraft.Figure 5 shows the raw reflectivity field, where a clearly defined hook echo can be distinguished in different radar volume scans (upper panels).The mesocyclone of the supercell is demonstrated with the radial velocity data (lower panels in Figure 5), where the typical couplet of velocity is depicted.In this case, the raw velocity data are overly thresholded; the filtering of the raw velocity data during the operational processing results in a loss of practically half of the mesocyclone in all the volumes (white pixels in dashed circles).As explained in Section 3.3, optimal tuning of the filters within SAMURAI allows us to fill data gaps during the analysis and therefore, to minimize the effect of the previous stringent operational filtering problem.The retrieved reflectivity and horizontal wind field at 1 km height is depicted in Figure 6a, where the southeasterly steering flow (Figure 3) is shown by the thick black arrow.A large area of weaker reflectivity values surrounds the main supercell, indicating that the supercell was embedded in the precipitation band, and the tornado was mostly rain-wrapped.In addition, Figure 6a depicts a region at the southeast corner where the wind is slightly more intense and converging within the rear part of the supercell, which is consistent with the storm inflow feeding the main updraft region (red shaded area in Figure 6c, indicating positive vertical velocity).The SAMURAI retrievals also resolve the closed circulation in the storm-relative horizontal winds (black arrows in Figure 6b), which depicts the mesocyclone at mid-levels (2.5 km AGL), and the area with a couplet of positive and negative vertical velocity (red and blue in Figure 6c).This vertical velocity couplet depicts the main updraft (red shaded) and the corresponding forward-flank downdraft (FFD; blue shaded) of this rain-wrapped heavy precipitation (HP) supercell [64][65][66].However, Figure 6c also depicts a complicated vertical vorticity retrieval, in 10 −5 s −1 , where the main couplet is displaced from the updraft-downdraft area, and therefore, we cannot completely relate the mesocyclone with the main updraft.Figure 7 shows the same retrieved SAMURAI fields for the supercell at 00:36 and 00:54 UTC, coinciding with the earlier stage entering the valley and when it produced the tornado damage on the ground, respectively.In both stages (Figure 7a,b), the closed circulation can be seen, indicating the mesocyclone and coinciding with the hook echo from Figure 5.However, again, we can see the decoupling of the mesocyclone from the main updraft, nearly collocated with the main downdraft, (Figure 7c,d), especially at 00:54 UTC when the tornado is on the ground and when the vertical velocity couplet is more intense (Figure 7d).This could be caused by the strong environmental shear strongly tilting the storm and decoupling the mesocyclone from the main updraft, especially at later stages, and this could be a realistic situation given the complexity of the case.However, it is also possible that this could be an effect of poor accuracy of the vorticity retrieval in this case, given how prone this field is to even qualitative error in any previous single-Doppler retrieval method, which could be misplacing the vorticity couplet.
Remote Sens. 2019, 11, x FOR PEER REVIEW 12 of 26 Figure 7 shows the same retrieved SAMURAI fields for the supercell at 00:36 and 00:54 UTC, coinciding with the earlier stage entering the valley and when it produced the tornado damage on the ground, respectively.In both stages (Figure 7a,b), the closed circulation can be seen, indicating the mesocyclone and coinciding with the hook echo from Figure 5.However, again, we can see the decoupling of the mesocyclone from the main updraft, nearly collocated with the main downdraft, (Figure 7c,d), especially at 00:54 UTC when the tornado is on the ground and when the vertical velocity couplet is more intense (Figure 7d).This could be caused by the strong environmental shear strongly tilting the storm and decoupling the mesocyclone from the main updraft, especially at later stages, and this could be a realistic situation given the complexity of the case.However, it is also possible that this could be an effect of poor accuracy of the vorticity retrieval in this case, given how prone this field is to even qualitative error in any previous single-Doppler retrieval method, which could be misplacing the vorticity couplet.Figure 8 shows the storm-relative horizontal (arrows) and vertical (blue contour lines) velocity at 00:36 and 00:54 UTC at 1.5 km and 2.5 km AGL over topography.It can be clearly seen how, at the earlier supercell stages, the wind flow is intense and channeled up the valley at 1.5 km (Figure 8a), and the mesocyclone, clearly seen at 2.5 km (Figure 8c) is also noticeable and located over the same region at lower levels (1.5 km, Figure 8a).Also, there's a strong correlation between the location of the low-level convergence area with the lower point of the valley (floor).The decoupling of the Figure 8 shows the storm-relative horizontal (arrows) and vertical (blue contour lines) velocity at 00:36 and 00:54 UTC at 1.5 km and 2.5 km AGL over topography.It can be clearly seen how, at the earlier supercell stages, the wind flow is intense and channeled up the valley at 1.5 km (Figure 8a), and the mesocyclone, clearly seen at 2.5 km (Figure 8c) is also noticeable and located over the same region at lower levels (1.5 km, Figure 8a).Also, there's a strong correlation between the location of the low-level convergence area with the lower point of the valley (floor).The decoupling of the mesocyclone with the main updraft is more evident at later stages, as seen in Figure 8b,d, where it is seen that the main updraft and downdraft are collocated with the windward and leeward side of the mountain (shaded in brown), which may be the cause of the enhancement of the updraft-downdraft pair at later stages as seen in Figure 7b,d.In this sense, the coupling of the strong environmental shear and the channeling up the valley, may have enhanced the FFD (converging winds towards the floor of the valley), and therefore caused stretching of the updraft and created the vertical vorticity at low-levels helping in the tornado formation, as well as playing a major role in the orientation and propagation of the supercell.Therefore, the complex supercell features retrieved with SAMURAI demonstrate that the technique was able to retrieve the storm-relative winds from heavily-filtered radar data, and that it is an effective tool to use with the Catalan operational C-band radar network.
Remote Sens. 2019, 11, x FOR PEER REVIEW 13 of 26 mesocyclone with the main updraft is more evident at later stages, as seen in Figure 8b,d, where it is seen that the main updraft and downdraft are collocated with the windward and leeward side of the mountain (shaded in brown), which may be the cause of the enhancement of the updraft-downdraft pair at later stages as seen in Figure 7b,d.In this sense, the coupling of the strong environmental shear and the channeling up the valley, may have enhanced the FFD (converging winds towards the floor of the valley), and therefore caused stretching of the updraft and created the vertical vorticity at lowlevels helping in the tornado formation, as well as playing a major role in the orientation and propagation of the supercell.Therefore, the complex supercell features retrieved with SAMURAI demonstrate that the technique was able to retrieve the storm-relative winds from heavily-filtered radar data, and that it is an effective tool to use with the Catalan operational C-band radar network.

Dynamical Influence of Topography
On 10 July 2013, a large splitting severe storm caused a substantial amount of losses of corn crops and fruit tree fields due to large hailstones (3-4 cm) and wind gusts in the Noguera county, located in the western Plains of Catalonia.The storm analyzed in [43] from a radar-based tracking point of view split into two new major cells, almost like mirror images (Figure 9a-c).In this splitting case, the raw radial velocity was not affected by such strong operational filtering (no missing bins in Figure 9d-f); however, the interaction of the storm with the local topography made this event challenging to monitor and nowcast.After the main split, the two resulting daughter storms (marked as numbers 2 and 3 in Figures 9 and 10) propagated completely in opposite directions (see arrows indicating LM

Dynamical Influence of Topography
On 10 July 2013, a large splitting severe storm caused a substantial amount of losses of corn crops and fruit tree fields due to large hailstones (3-4 cm) and wind gusts in the Noguera county, located in the western Plains of Catalonia.The storm analyzed in [43] from a radar-based tracking point of view split into two new major cells, almost like mirror images (Figure 9a-c).In this splitting case, the raw radial velocity was not affected by such strong operational filtering (no missing bins in Figure 9d-f); however, the interaction of the storm with the local topography made this event challenging to monitor and nowcast.After the main split, the two resulting daughter storms (marked as numbers 2 and 3 in Figures 9 and 10) propagated completely in opposite directions (see arrows indicating LM and RM in Figure 10).The split resulted in a long-lasting LM, which propagated for more than two hours across the topography towards northeast.Finally, after the main splitting, the RM merged with a new cell travelling from the northwest (cell 4 in Figure 10), resulting in an upscale growth of a multicell system that travelled southeastward (cell 5 in Figure 10).The convective line generated from that multicell system resulted in severe wind gusts, indicated by red stars near cell 5 in Figure 10.Red stars indicate severe weather in the form of hail or damaging wind gusts that the sparse spotters could report in an isolated area.It is important to notice that the wind gusts could be felt several kilometers away from the storm centroid (southeast red star in Figure 10), due to the size of the storm and the strong outflow from the first splitting storm.
One of the primary ingredients needed for supercells to split is strong deep-layer environmental shear.A common synoptic-scale flow that favors such conditions is a strong jet at 300 hPa with the exit region over the storm area, enhancing the instability and lifting conditions, along with a tilted and strong low-level jet (850 hPa) that may enhance the shear (see, for instance, [67][68][69]).However, on 10 July 2013, the left exit region of the jet at 300 hPa was displaced from the convection in Catalonia and was located over the eastern Mediterranean Sea due to a cut-off low at that same level north of Catalonia.This created a strong northwesterly upper-level right jet entrance region located over the northern coast of Spain, placing Catalonia at the left entrance of the jet streak and therefore, not in the maximum vorticity advection area for optimal lifting [70].and RM in Figure 10).The split resulted in a long-lasting LM, which propagated for more than two hours across the topography towards northeast.Finally, after the main splitting, the RM merged with a new cell travelling from the northwest (cell 4 in Figure 10), resulting in an upscale growth of a multicell system that travelled southeastward (cell 5 in Figure 10).The convective line generated from that multicell system resulted in severe wind gusts, indicated by red stars near cell 5 in Figure 10.
Red stars indicate severe weather in the form of hail or damaging wind gusts that the sparse spotters could report in an isolated area.It is important to notice that the wind gusts could be felt several kilometers away from the storm centroid (southeast red star in Figure 10), due to the size of the storm and the strong outflow from the first splitting storm.One of the primary ingredients needed for supercells to split is strong deep-layer environmental shear.A common synoptic-scale flow that favors such conditions is a strong jet at 300 hPa with the exit region over the storm area, enhancing the instability and lifting conditions, along with a tilted and strong low-level jet (850 hPa) that may enhance the shear (see, for instance, [67][68][69]).However, on 10 July 2013, the left exit region of the jet at 300 hPa was displaced from the convection in Catalonia and was located over the eastern Mediterranean Sea due to a cut-off low at that same level north of Catalonia.This created a strong northwesterly upper-level right jet entrance region located over the northern coast of Spain, placing Catalonia at the left entrance of the jet streak and therefore, not in the maximum vorticity advection area for optimal lifting [70].
The forecasted WRF sounding used in the analysis for this case is shown in Figure 11.Lleida was the closest WRF sounding location, 36 km southwest of the storm location (Figure 11b), valid at 12 UTC for a WRF run initialized at 00 UTC.Although not strictly resulting from the jet streak, the northeasterly flow over Catalonia was significantly strong at high levels (Figure 11), which helped build moderate deep layer shear (22 kt), since the wind in the boundary layer was negligible (Figure 11a).The deep layer shear is more conducive of multicell systems rather than supercells [61,65], which explains the behavior of the storm at the beginning.In addition, it can be seen that, in the case of splitting, the Bunkers storm motion vectors [62] depicted an enhanced LM towards the northeast, coinciding with the final motion after the split.Nevertheless, the 12 UTC WRF sounding does not fully explain the rotation of the resulting cells, thus the local topography and its orientation relative to the storm could have played an important role.The forecasted WRF sounding used in the analysis for this case is shown in Figure 11.Lleida was the closest WRF sounding location, 36 km southwest of the storm location (Figure 11b), valid at 12 UTC for a WRF run initialized at 00 UTC.Although not strictly resulting from the jet streak, the northeasterly flow over Catalonia was significantly strong at high levels (Figure 11), which helped build moderate deep layer shear (22 kt), since the wind in the boundary layer was negligible (Figure 11a).The deep layer shear is more conducive of multicell systems rather than supercells [61,65], which explains the behavior of the storm at the beginning.In addition, it can be seen that, in the case of splitting, the Bunkers storm motion vectors [62] depicted an enhanced LM towards the northeast, coinciding with the final motion after the split.Nevertheless, the 12 UTC WRF sounding does not fully explain the rotation of the resulting cells, thus the local topography and its orientation relative to the storm could have played an important role.
The retrieved reflectivity and the horizontal wind field (upper and mid-row), as well as the vertical velocity and vertical absolute vorticity (lower row), before and after the splitting, are shown in Figure 12.First, it is worth mentioning that the splitting is not clearly visible until 16:06 UTC (Figure 12b), and even at that time, the shape does not coincide with the usual splitting pattern.During the time steps before splitting (Figure 12a), the entire pattern is more similar to a multicell system, since during almost one hour, different small cells grow and dissipate near each other (see, for instance Figure 12b,d,e).Figure 12a,b show a strong shear in the storm area, with southeasterly flow at low levels (Figure 12a,b), and a strong northwesterly flow at upper levels (Figure 12d,e), consistent with the synoptic situation at 300 hPa commented on previously.Figure 12d,e show that at Z = 4 km, there is a formation of convergence lines of winds (box and dashed lines labeled CL), indicating another region where new cells are formed.In a heterogeneous multicell environment, that configuration would help the formation of the new cells upstream of the low-level flow, resulting from the lifting mechanism of the gust front of the preceding cells (downshear flank of the outflow).This is visible in Figure 12h marked with a dashed semicircle, where it can be seen how different updrafts appear upstream of the low-level wind, although not long-lasting.Despite that clear configuration, how the general movement of the cells is towards the southwest is also seen (note, the change in the y-component of the storm location between Figure 12d,e).The retrieved reflectivity and the horizontal wind field (upper and mid-row), as well as the vertical velocity and vertical absolute vorticity (lower row), before and after the splitting, are shown in Figure 12.First, it is worth mentioning that the splitting is not clearly visible until 16:06 UTC (Figure 12b), and even at that time, the shape does not coincide with the usual splitting pattern.During the time steps before splitting (Figure 12a), the entire pattern is more similar to a multicell system, since during almost one hour, different small cells grow and dissipate near each other (see, for instance Figure 12b,d,e).Figure 12a,b show a strong shear in the storm area, with southeasterly flow at low levels (Figure 12a,b), and a strong northwesterly flow at upper levels (Figure 12d,e), consistent with the synoptic situation at 300 hPa commented on previously.Figure 12d,e show that at Z = 4 km, there is a formation of convergence lines of winds (box and dashed lines labeled CL), indicating another region where new cells are formed.In a heterogeneous multicell environment, that configuration would help the formation of the new cells upstream of the low-level flow, resulting from the lifting mechanism of the gust front of the preceding cells (downshear flank of the outflow).This is visible in Figure 12h marked with a dashed semicircle, where it can be seen how different updrafts appear upstream of the low-level wind, although not long-lasting.Despite that clear configuration, how the general movement of the cells is towards the southwest is also seen (note, the change in the y-component of the storm location between Figure 12d,e).The retrieved horizontal wind field and reflectivity over topography is shown in Figure 13 to better illustrate the impact of topography on the storm behavior and propagation.First, at the earliest stages of the initial storm, there is a divergence of the wind (red arrows in Figure 13d) and a convergence zone at mid-levels (blue arrows in Figure 13d), coinciding with the orientation of the main ridges in the valley (notice the convergence areas parallel to the southwest lower ridge, SWR, and the divergence of the north coinciding with the entrance of the wind between the two northwest higher ridges, NWR).The role of the northwest ridges (NWR) is more visible in later stages of the storm.Figure 13b shows again the strong south-southeasterly low-level winds impinging against the northwest ridges (framed in black in Figure 13b,c), coincident with the gust front from Figure 12h,i, which demonstrate the role of the northern ridges (NWR) in the triggering of new cells in the downshear flank, and probably, in the enhancement of the preexisting cells.The coupling of those new cells and the strong low-level flow towards the northwest led to faster propagation of the entire system towards that direction and enhanced the strength of cores in that region.In addition, Figure 13e shows the storm starting to split at upper levels (notice that the two reflectivity cores are more defined in Figure 13e than Figure 13b), placing the two main cores (LM and RM in Figure 13e) on either side of the southwest ridge (SWR).The start of the mesocyclone formation at mid-levels is shown with an orange circle in Figure 13e, coincident with a convergence area at 1.5 km (orange circle in Figure 13b), where the winds seem to be channeled, enhancing the main updraft of the RM supercell.The orientation of the reflectivity field in relation to the valley and the main ridges is more clear at later stages, as seen in Figure 13c,f.It can be seen how the two main supercells are almost identical to the shape of the valley (Figure 13c), taking an elongated shape towards the northwest.Finally, it is seen how at mid-levels, the difluence of the wind is even stronger over the NWR ridge (red frame in Figure 13f).The NWR ridge separates the entire system, creating two major flanks, probably forcing them to move in opposite directions.Additionally, the winds channeling through the tallest northeast ridges (NER) (blue arrows in Figure 13f) create different converging regions at mid-levels, which could have enhanced the LM updraft, being one of the causes for its intensification and its long-lasting motion.The retrieved horizontal wind field and reflectivity over topography is shown in Figure 13 to better illustrate the impact of topography on the storm behavior and propagation.First, at the earliest stages of the initial storm, there is a divergence of the wind (red arrows in Figure 13d) and a convergence zone at mid-levels (blue arrows in Figure 13d), coinciding with the orientation of the main ridges in the valley (notice the convergence areas parallel to the southwest lower ridge, SWR, and the divergence of the north coinciding with the entrance of the wind between the two northwest higher ridges, NWR).The role of the northwest ridges (NWR) is more visible in later stages of the storm.Figure 13b shows again the strong south-southeasterly low-level winds impinging against the northwest ridges (framed in black in Figure 13b,c), coincident with the gust front from Figure 12h and i, which demonstrate the role of the northern ridges (NWR) in the triggering of new cells in the downshear flank, and probably, in the enhancement of the preexisting cells.The coupling of those new cells and the strong low-level flow towards the northwest led to faster propagation of the entire system towards that direction and enhanced the strength of cores in that region.In addition, Figure 13e shows the storm starting to split at upper levels (notice that the two reflectivity cores are more   Figure 14 shows horizontal reflectivity and wind at 4 km before and after the split of the storm, as well as the vertical cross sections along the different lines (A-B, C-D, and E-F).The cross section in Figure 14a, along the southwest-northeast direction, depicts the formation of the new cells at mid-levels ~6 km AGL (6 km to the east of point A in Figure 14b), where a new cell core can be seen with a strong lifting in the forward flank (updraft indicated with solid contour lines, Figure 14b).This new cell corresponds to point (−20, 27) in Figure 14a.Although this configuration seems to be predominant in the initial time steps, causing the new cells to converge in the southwest region and creating a major flank, it follows a similar pattern in the opposite direction (northeast), where another convergence line is created and a major cell starts to form.This is seen with the inverted-U shape in the reflectivity contours in Figure 13e.Figure 14e shows the tilting of the entire system towards the southwest (along C-D axis in Figure 14e), which also depicts the origin of the main two updrafts and a possible starting RM rotation (vorticity framed in black in Figure 14f).The splitting process then occurs between 16:00 UTC and 16:12 UTC (not shown), where the two main cores organize the creating of the two-vortex couplet visible at 16:24 UTC (Figure 14l).The result is two supercells (Figure 14j,k) moving in opposite directions with cyclonic and anticyclonic circulations.Figure 14l, which shows the storm cross section along the E-F axis, clearly depicts the two main cells, completely separated and tilted towards the direction of motion.Two couplets of updrafts-downdraft from mid-upper levels can also be seen.This could be related to the main inflow direction and the rear-flank downdraft.Additionally, Figure 14l shows two main couplets in this case of vertical vorticity.It is seen how the maxima (900 × 10 −5 s −1 ) and minima (−300 × 10 −5 s −1 ) values at 16:24 UTC are located around 4 km, depicting the resulting LM and RM supercells, as expected from the splitting process.

Discussion and Conclusions
In the present study, an open-source variational multi-Doppler software, SAMURAI, was applied in Catalonia (northeast of the Iberian Peninsula).The area of study presents complex topography and is prone to severe weather events and heavy rainfall.This work constitutes a pioneering study in a southern European country, thanks to the geographical configuration of the

Discussion and Conclusions
In the present study, an open-source variational multi-Doppler software, SAMURAI, was applied in Catalonia (northeast of the Iberian Peninsula).The area of study presents complex topography and is prone to severe weather events and heavy rainfall.This work constitutes a pioneering study in a southern European country, thanks to the geographical configuration of the XRAD, the radar network of the SMC.It is the first time that SAMURAI has been used for analyzing severe weather at a thunderstorm scale and in complex terrain using SMC operational radar data.This paper illustrates the capabilities of the software for solving different situations, in which the anomalies of the original data allow for characterization of the internal thunderstorm dynamics due to aliasing and aggressive operational velocity filters removing part of the data, especially in high azimuthal shear regions.These data issues are common in most of the operational European weather radar networks.
With the analysis of the tornadic thunderstorm on 7 January 2018, the present work demonstrates that SAMURAI constitutes a good tool for enabling storm analysis in complex terrain: a tornado-producing supercell channeling through a valley.The SAMURAI analyses are promising because they seem to reasonably retrieve the horizontal and vertical (especially at mid-levels) storm-related wind field even though i) there is radar beam blockage and the lowest radar beam scans above the tornado base, ii) the radars are operational with a fixed radar scan mode that may result in underrepresentation of the storm winds at upper levels, and iii) aggressive pre-filtering of raw radial velocity data, resulting in a loss of half of the mesocyclone.By using a spline-based variational method that includes low-pass filtering, we have been able to retrieve the mid-level mesocyclone of the tornadic thunderstorm, as well as the expected updraft-downdraft couplet associated with the circulation.
Although large qualitative errors are typical in variational single-Doppler analysis, in the case of multi-Doppler retrievals, variational algorithms such as SAMURAI have been demonstrated to minimize such errors [21].Considering that our analyses are made with two radars, we can consider the solutions found in this research to be reliable.However, we need to bear in mind the possible uncertainties that can arise in the retrievals.For instance, the smoothing and constraints of the data in the variational analysis, the signal attenuation, or the constant advection flow use for the storm motion can impact the wind retrievals (e.g., [71,72]).These uncertainties are larger at upper levels, where the radar data are less dense, and it has been recently demonstrated that the radar scanning strategy can underestimate strong updrafts at low levels and overestimate them in upper levels [73].While the way SAMURAI handles mass continuity via the local application of the constraint in derivative form rather than via integration of the equation itself helps to reduce these uncertainties, it is difficult to overcome the Doppler geometry limitations [74] completely, and the imperfect boundary conditions and missing low-level divergence will still affect the derived winds, especially when working with complex fine-scale flows.The need to improve how SAMURAI handles the winds near the surface still remains, particularly when the storm is over topography.However, this work demonstrates that SAMURAI can be applied to retrieve physically reasonable mid-level circulations and provide information about key features of organized thunderstorms, such as mesocyclones, updrafts, and vorticity couplets, even while using an operational radar network with a fixed scanning strategy.The software is currently being improved by specifying directly the topography in the analysis, which will enhance the retrieved near-surface winds in complex terrain in future analyses.
The analysis of the anomalous (weak convective environment) splitting thunderstorm case on 10 July 2013 has shown that the technique allows for the identification of the features and internal core interactions that lead to the separation of the two main cells.In the analyzed case, we have detected that, contrary to most of the splitting storms reported in the literature that present strong shear in the environmental winds and a primary single cell, the final splitting process here reported was a response of a major system (multicell) interacting with its surrounding conditions, especially the topography.This was the result of an enhancement and organization of different cores in the rear flank of a multicell system, due to the proximity to a steep orographic ridge, which also enhanced the storm-related mid-upper-level winds.The splitting process occurred after the organization of the main two cells, which started at high levels, and no initial vertical vorticity couplet was depicted before the split.SAMURAI retrievals provided information about the vertical structure, especially at mid-levels, and the associated wind field of the thunderstorms during the entire life cycle.
Although the current SAMURAI version has been optimized to run five times faster than the first LROSE Blaze release (used in this research), the entire process for retrieving the 3D winds for a radar volume takes more than six minutes (operational time at the SMC).This means that the analysis cannot yet be applied to keep up with real-time.For this reason, the retrievals are purely for research purposes and post-event analysis, which is still beneficial for the Catalan forecasters, who can work with better knowledge in similar situations affecting the area.However, the LROSE team is working towards improving the Fast Reorder and CEDRIC Technique in LROSE (FRACTL) algorithm.This is an optimized version of the original Reorder gridding [49] and Custom Editing and Display of Reduced Information in Cartesian space (CEDRIC, [75]) classic dual-Doppler retrieval, which could be used operationally, since the analysis is done much faster.In fact, our first FRACTL investigations were conducted on the NCAR-Earth Observing Laboratory (EOL) computers for research purposes and took only three minutes to retrieve an entire volume.Unfortunately, the preliminary FRACTL investigations did not produce realistic results, possibly due to the use of a classic retrieval instead of a variational solution with the Catalan radar data and the complex topography, and it is not yet clear that this could improve operational nowcasting.Furthermore, we expect that the same work at the SMC would require a longer time due to the SMC computational skills (regional weather service).In this sense, the SMC is working towards improving the aggressive filtering of the data and on recording the original unfiltered radar moments, which, along with the improvements in the FRACTL and SAMURAI retrievals, may lead Catalonia to the possibility of having operational multi-Doppler retrievals in the future.The use of this methodology in the cases presented, with complicated internal dynamics, results in an improvement in the knowledge of the nature of these anomalous thunderstorms which are frequent in Catalonia.Furthermore, with this analysis, it is demonstrated that in order to understand how severe thunderstorms evolve and behave during the entire life cycle, it is necessary to put effort into retrieving the wind field at all levels, not just to focus on what happens near the surface.This is usually achieved with idealized numerical simulations at high resolution or with multiple-Doppler analyses in field campaigns, with X-band radars providing better spatial resolution than operational radar networks.However, for the first time in Catalonia (and also, the Iberian Peninsula), complete storm-relative wind fields have been obtained purely from operational networks, and in particular, with an operational C-band radar network.This has provided a three-dimensional view of the wind field inside a thunderstorm, constituting an important step forward in the knowledge of topographic influence on western Mediterranean thunderstorms We plan to consider other events in the coastal areas, where the dynamics can be affected by other factors, such as sea breezes, mesoscale convergence lines, or maritime boundaries and not only by the internal dynamic of the thunderstorms or interaction with other convective systems and the complex topography.The Mediterranean coast has become a hotspot for heavy rains and severe weather in recent years, and future scenarios demonstrate that the frequency of such events will increase with global warming (see, for instance, [76,77]).There is a current research focus in the Mediterranean basins, especially in densely populated areas, that require concrete action plans and improved early warning systems for future adaptation [5].The tools presented here and the configuration of the XRAD will allow the Catalan community to further study the convection affecting coastal areas, which may also benefit other risk areas such as Italy or the coast of Greece, where there is not a good radar network configuration to retrieve multi-Doppler analyses.

Figure 1 .
Figure 1.Geographical distribution of severe storm reports in Catalonia.(a) Density of severe weather events per county (events•km −2 ).(b) Density of population per county (inhabitants•km −2 ).Counties are indicated in grey lines.Legend categories are based on quantiles.

Figure 1 .
Figure 1.Geographical distribution of severe storm reports in Catalonia.(a) Density of severe weather events per county (events•km −2 ).(b) Density of population per county (inhabitants•km −2 ).Counties are indicated in grey lines.Legend categories are based on quantiles.

Figure 2 .
Figure 2. The Meteorological Service of Catalonia radar network, XRAD.(a) Spatial coverage (shaded in grey) and topography blockage (areas in black) for 0.6 elevation angle.All radars have a maximum range of 130 km, except for CDV (central one) which is 150 km.(b) Dual-Doppler pairs for the XRAD: CDV-PBE and CDV-LMI.Dotted circles show the dual lobes for a 30 dual-Doppler beam-crossing angle.Topography is shaded in terrain colors and is indicated in meters.The main topographic features are labeled.

Figure 2 .
Figure 2. The Meteorological Service of Catalonia radar network, XRAD.(a) Spatial coverage (shaded in grey) and topography blockage (areas in black) for 0.6 • elevation angle.All radars have a maximum range of 130 km, except for CDV (central one) which is 150 km.(b) Dual-Doppler pairs for the XRAD: CDV-PBE and CDV-LMI.Dotted circles show the dual lobes for a 30 • dual-Doppler beam-crossing angle.Topography is shaded in terrain colors and is indicated in meters.The main topographic features are labeled.

Figure 3 .
Figure 3. (a) Observed atmospheric profile from the Barcelona sounding on 7 January 2018 at 00 UTC.The hodograph shows wind velocity and direction up to 6 km height, in gradient blue, with dark blue starting from the surface.Potential left-(LM) and right-moving (RM) Bunkers supercell motions [62] are depicted by black dots.The black dot at the surface represents the Lifting Condensation Level (LCL).Sounding retrieved with MetPy Python's module [63].(b) Barcelona sounding location with respect to the tornado-producing storm track.

Figure 3 .
Figure 3. (a) Observed atmospheric profile from the Barcelona sounding on 7 January 2018 at 00 UTC.The hodograph shows wind velocity and direction up to 6 km height, in gradient blue, with dark blue starting from the surface.Potential left-(LM) and right-moving (RM) Bunkers supercell motions [62] are depicted by black dots.The black dot at the surface represents the Lifting Condensation Level (LCL).Sounding retrieved with MetPy Python's module [63].(b) Barcelona sounding location with respect to the tornado-producing storm track.

Figure 4 .
Figure 4. Track of the nocturnal winter tornado-producing storm on Central Catalonia on January 7 2018 over topography.Part of the 30 CDV-PBE pair dual lobes is depicted in the dotted lines.Dots indicate supercell centroid times labeled in UTC. Green dots indicate centroids where the mesocyclone was depicted in the Doppler velocity field as a radial velocity couplet.Red stars indicate confirmed tornado damage on the ground (post-event analysis).The black arrow indicates the steering flow at 3 km.

Figure 4 .
Figure 4. Track of the nocturnal winter tornado-producing storm on Central Catalonia on January 7 2018 over topography.Part of the 30 • CDV-PBE pair dual lobes is depicted in the dotted lines.Dots indicate supercell centroid times labeled in UTC. Green dots indicate centroids where the mesocyclone was depicted in the Doppler velocity field as a radial velocity couplet.Red stars indicate confirmed tornado damage on the ground (post-event analysis).The black arrow indicates the steering flow at 3 km.

Figure 6 .
Figure 6.Retrieved SAMURAI fields for the supercell case on 7 January 2018 at 00:42 UTC.The thick black arrow in upper left indicates the steering flow at 3 km, and the black frame indicates the storm inflow.Reflectivity (dBZ) and storm-relative horizontal wind field (m•s −1 ) with black arrows at (a) Z = 1.0 km and (b) Z = 2.5 km.(c) Vertical velocity in color scale (w, in m•s −1 ) and vertical vorticity (ζ, in 10 −5 •s −1 ), represented with black contours (dashed and solid for negative and positive values, respectively).Contours are labeled every 200 × 10 −5 •s −1 (0 value not shown).The 35 dBZ contour is depicted with the thick black line.

Figure 6 .
Figure 6.Retrieved SAMURAI fields for the supercell case on 7 January 2018 at 00:42 UTC.The thick black arrow in upper left indicates the steering flow at 3 km, and the black frame indicates the storm inflow.Reflectivity (dBZ) and storm-relative horizontal wind field (m•s −1 ) with black arrows at (a) Z = 1.0 km and (b) Z = 2.5 km.(c) Vertical velocity in color scale (w, in m•s −1 ) and vertical vorticity (ζ, in 10 −5 •s −1 ), represented with black contours (dashed and solid for negative and positive values, respectively).Contours are labeled every 200 × 10 −5 •s −1 (0 value not shown).The 35 dBZ contour is depicted with the thick black line.

Figure 6 .
Figure 6.Retrieved SAMURAI fields for the supercell case on 7 January 2018 at 00:42 UTC.The thick black arrow in upper left indicates the steering flow at 3 km, and the black frame indicates the storm inflow.Reflectivity (dBZ) and storm-relative horizontal wind field (m•s −1 ) with black arrows at (a) Z = 1.0 km and (b) Z = 2.5 km.(c) Vertical velocity in color scale (w, in m•s −1 ) and vertical vorticity (ζ, in 10 −5 •s −1 ), represented with black contours (dashed and solid for negative and positive values, respectively).Contours are labeled every 200 × 10 −5 •s −1 (0 value not shown).The 35 dBZ contour is depicted with the thick black line.

Figure 8 .
Figure 8. Retrieved SAMURAI fields for the supercell case on 7 January 2018 at 00:36 and 00:54 UTC.Storm-relative horizontal wind field (m•s −1 , with black arrows) and vertical velocity (w, every 3 m•s −1 , with blue contours, dashed and solid for negative and positive values, respectively) at (a,b) Z = 1.5 km and (c,d) Z = 2.5 km.

Figure 8 .
Figure 8. Retrieved SAMURAI fields for the supercell case on 7 January 2018 at 00:36 and 00:54 UTC.Storm-relative horizontal wind field (m•s −1 , with black arrows) and vertical velocity (w, every 3 m•s −1 , with blue contours, dashed and solid for negative and positive values, respectively) at (a,b) Z = 1.5 km and (c,d) Z = 2.5 km.

Figure 10 .
Figure 10.Track of the splitting storm 10 July 2013 over topography.Part of the 30° CDV-LMI pair dual-lobes is depicted in dotted lines.Red stars indicate Severe Weather (SW) occurrences (wind gusts and hail).Red numbers indicate the different cells that took place in the system identified by [43], being 1 (primary cell), 2 (left mover), and 3 (right mover), the ones of interest taking part in the splitting process.Red squares indicate the starting point of the left (cell 2) and right (cell 3) movers, and red arrows indicate their motion.The black arrow indicates the mean NW steering flow at 3 km (direction of motion for the majority of the cells during that day).Dots indicate the different cell centroid times labeled in UTC.

Figure 10 .
Figure 10.Track of the splitting storm 10 July 2013 over topography.Part of the 30 • CDV-LMI pair dual-lobes is depicted in dotted lines.Red stars indicate Severe Weather (SW) occurrences (wind gusts and hail).Red numbers indicate the different cells that took place in the system identified by [43], being 1 (primary cell), 2 (left mover), and 3 (right mover), the ones of interest taking part in the splitting process.Red squares indicate the starting point of the left (cell 2) and right (cell 3) movers, and red arrows indicate their motion.The black arrow indicates the mean NW steering flow at 3 km (direction of motion for the majority of the cells during that day).Dots indicate the different cell centroid times labeled in UTC.

26 Figure 11 .
Figure 11.(a) Forecasted WRF sounding for Lleida on 10 July 2013, initialized at 00 UTC and valid at 12 UTC.The hodograph shows the wind velocity and direction up to 6 km height, in gradient blue, with dark blue starting from the surface.Potential left-(LM) and right-moving (RM) Bunkers supercell motions [62] are depicted with black dots.The black dot at ~750 hPa represents the Lifting Condensation Level (LCL).Sounding retrieved with MetPy Python's module [63].(b) Lleida WRF sounding location with respect to storm track.

Figure 11 .
Figure 11.(a) Forecasted WRF sounding for Lleida on 10 July 2013, initialized at 00 UTC and valid at 12 UTC.The hodograph shows the wind velocity and direction up to 6 km height, in gradient blue, with dark blue starting from the surface.Potential left-(LM) and right-moving (RM) Bunkers supercell motions [62] are depicted with black dots.The black dot at ~750 hPa represents the Lifting Condensation Level (LCL).Sounding retrieved with MetPy Python's module [63].(b) Lleida WRF sounding location with respect to storm track.

26 Figure 12 .
Figure 12.SAMURAI retrievals for the splitting case on 10 July 2013 at Z = 0.5 km (upper row), and Z = 4 km (mid and lower rows).(a-f) Reflectivity (dBZ) and storm-relative horizontal wind field (m•s −1 ) (with black arrows).(g-i) Vertical velocity (w, in m•s −1 ) and vertical vorticity (ζ, in 10 −5 •s −1 ), represented with black contours (dashed and solid for negative and positive values, respectively).Contours are labeled every 200 × 10 −5 •s −1 .Boxes indicate areas of convergence winds, and dashed lines indicate the axis of the convergence line formed by those winds.CL and GF labels indicate convergence line and gust front, respectively.Black arrows indicate the two final vortex couplets.

Figure 12 .
Figure 12.SAMURAI retrievals for the splitting case on 10 July 2013 at Z = 0.5 km (upper row), and Z = 4 km (mid and lower rows).(a-f) Reflectivity (dBZ) and storm-relative horizontal wind field (m•s −1 ) (with black arrows).(g-i) Vertical velocity (w, in m•s −1 ) and vertical vorticity (ζ, in 10 −5 •s −1 ), represented with black contours (dashed and solid for negative and positive values, respectively).Contours are labeled every 200 × 10 −5 •s −1 .Boxes indicate areas of convergence winds, and dashed lines indicate the axis of the convergence line formed by those winds.CL and GF labels indicate convergence line and gust front, respectively.Black arrows indicate the two final vortex couplets.

Figure 13 .
Figure 13.SAMURAI retrievals over topography for the splitting case on 10 July 2013 at Z = 1.5 km (upper row), and Z = 4 km (lower row).Reflectivity (dBZ) in contour lines (labeled every 5 dBZ) and storm-relative horizontal wind field (m•s −1 ) (with black arrows) at (a-d) 15:48 UTC, (b-e) 16:12 UTC, and (c-f) 16:24 UTC.NWR, NER, and SWR indicate the northwest, northeast, and southwest ridges, respectively.Black and red frames indicate regions with winds impinging against the NWR ridges, and the region with wind difluence at mid-levels, respectively.Thick blue and red arrows indicate areas with converging and diverging winds, respectively.Orange circles indicate the closed circulation for the left mover supercell.Finally, LM and RM indicate the left-and right-mover supercells.

Figure 13 .
Figure 13.SAMURAI retrievals over topography for the splitting case on 10 July 2013 at Z = 1.5 km (upper row), and Z = 4 km (lower row).Reflectivity (dBZ) in contour lines (labeled every 5 dBZ) and storm-relative horizontal wind field (m•s −1 ) (with black arrows) at (a-d) 15:48 UTC, (b-e) 16:12 UTC, and (c-f) 16:24 UTC.NWR, NER, and SWR indicate the northwest, northeast, and southwest ridges, respectively.Black and red frames indicate regions with winds impinging against the NWR ridges, and the region with wind difluence at mid-levels, respectively.Thick blue and red arrows indicate areas with converging and diverging winds, respectively.Orange circles indicate the closed circulation for the left mover supercell.Finally, LM and RM indicate the left-and right-mover supercells.

Table 1 .
Radar site name, abbreviation, height AGL (m), location, operational unambiguous maximum ranges (km) and operational elevation angles ( • ), total number of elevations, and gate spacing (m) of the XRAD.

Table 2 .
Dual-Dopplerri network coverage characteristics for the XRAD, considering a crossing-beam angle of 30 • .The two pairs used are shaded in grey.