A versatile scanning method for volumetric measurements of velocity and density fields

Understanding turbulence in a stratified environment requires a detailed picture of both the velocity field and the density field. Experimentally, this represents a significant measurement challenge, especially when full three-dimensional data is needed to accurately characterise the turbulent fields. This paper presents a new approach to obtaining such data through well-resolved, near-instantaneous volume-spanning measurements. This is accomplished using a simultaneous particle image velocimetry and planar-laser-induced fluorescence scanning setup that is highly versatile while allowing precise overlapping light sheets in a volume. Moreover, the approach detailed here allows for complex scan sequences over varying volumes without the need for mechanical adjustment. To demonstrate the technique, preliminary results from the buoyancy-driven exchange flow through an inclined duct connecting two reservoirs of different density fluids are shown.


Introduction
Flows of geophysical and environmental interest are frequently turbulent, exhibiting a range of scales and chaotic motions characterised by a large Reynolds number Re = UL ν 1.
Here, U is a typical velocity scale, L is a typical length scale and ν is the kinematic viscosity. Moreover, buoyancy forces (due to density differences) play an important role in these flows and are often characterised by a Richardson number Ri = N 2 L 2 U 2 , where N is the buoyancy frequency. Despite advances in experimental measurement techniques, capturing near-instantaneous, simultaneous volumetric density and velocity fields remains challenging. Nevertheless, having such information is desirable in the study of stratified turbulence due to the three-dimensional (3D) nature of the turbulent flow and the dynamical significance of the density field. Such data would also be helpful if coherent structures are to be extracted from the flow. Coherent structures have offered great insight into unstratified turbulent flow (e.g. wall-bounded turbulence (Jiménez 2018)) and, more recently, stratified flow (e.g. Lucas et al (2017)). As well as this, there is still a need for experimental measurements of realisable flows to complement direct numerical simulations (DNS). This is due to the complexity of modelling certain physical boundary conditions and domains, as well as the non-trivial dependence on the Schmidt number Sc = ν κ , where κ is the molecular diffusivity. As demonstrated by Zhou et al (2017), the behaviour of stratified flows depend on Sc. However, DNS of high-Sc stratified flows (Sc = 10 2 -10 3 in the salinity-stratified ocean) require sufficient resolution to not only capture the Kolmogorov scale, but also the Bachelor scale, respectively defined as η k = ν 3 1/4 and η b = η k where is the dissipation rate of turbulent kinetic energy. Currently, this resolution requirement makes high-Sc DNS prohibitive, even at the relatively low Re and that are easily achieved using salt-stratified (high-Sc) laboratory experiments. To this end, we present a new approach to measure simultaneous volumetric density and velocity fields at the high Re and Sc typical of laboratory flows. The layout of the remainder of this paper is as follows. In section 2 we describe some of the challenges associated with the current methodologies to measure velocity and density fields in stratified flows before describing our new approach in section 3. Our experimental setup, including the specific hardware used to perform the measurements, is detailed in section 4. Illustrative results are presented in section 5. Finally, discussions and conclusions are given in section 6.

Velocity measurements
Here we discuss some of the challenges associated with velocity measurement techniques. Two volumetric methods are compared: non-scanned and scanned. In non-scanned approaches the velocity field is calculated by illuminating and imaging (with multiple cameras) the whole volume of interest instantaneously. This differs from scanned techniques, where planar measurements calculated from thin light sheets are combined to build a volume. A full review of the techniques is beyond the scope of this paper and the reader is referred to Discetti and Coletti (2018) for more information.
A challenge when making optical measurements such as particle tracking velocimetry (PTV) or particle image velocimetry (PIV) in stratified flows is dealing with refractive index variations. The refractive index variations can be minimised by choosing two solutions that have different densities but the same refractive index (e.g. McDougall (1979), Dalziel et al (1999) and Verso et al (2017)), but inevitably small residual refractive index variations remain. Unfortunately, nonscanned volumetric approaches are more susceptible to any residual refractive index mismatch than scanning measurements, making them less suited to the study of stratified flows. This is due to non-scanned approaches relying on accurately triangulating the particles' 3D positions and calculating their 3D displacements at subsequent times. In comparison, Dalziel et al (2007) demonstrated that the planar PIV approach, used in the scanning measurements, is more robust to such errors as the pattern matching algorithm is less sensitive to absolute refractive index variations than to their rate of change.
As well as the refractive index issues, spatial resolution is also a limiting factor for non-scanned approaches as all of the particles within the measurement volume have to be distinguishable simultaneously in the camera images. This means that the yield of vectors in the total volume is limited, with either the spatial resolution or spatial extent being compromised. The resolution of a non-scanned system is typically expressed in particles per pixel (ppp). The latest algorithms can achieve 0.1 ppp, yielding one velocity vector per particle (Schanz and Schröder 2016). To compare this with a scanning system, consider a resolution of one velocity vector every 8 × 8 pixels, a four megapixel camera will give ∼ 6 × 10 4 vectors on the plane for PIV compared with ∼ 4 × 10 5 vectors for the 0.1 ppp non-scanned system. However, by building a volume of more than six planes, the scanning approach produces a higher volumetric vector yield. On the other hand, an advantage to the non-scanning approach is that the 3D field is effectively instantaneous: there is no time lag when imaging across the volume and the flow is adequately frozen by the pulse width of the light source.
Scanning methods have their own unique challenges. The fundamental issue when using the scanning approach is ensuring sufficient light sheet overlap between a pair of images to permit accurate velocity measurements, i.e. keeping sufficient particle overlap between image pairs. Previous work by Brücker (1997) overcame this issue with the novel use of mirrors on a rotating drum to obtain a series of discrete pairs of light sheets that were coincident in space. However, the use of such a mechanism relies on mechanical reconfiguration to alter the scanning volume and does not lend itself to more complicated scan sequences. More recent studies have opted to scan a light sheet using oscillating mirrors positioned using galvanometers (for accurate and fast positioning). For example, Krug et al (2014) produced a light sheet by passing a laser beam through a cylindrical lens before scanning the sheet using an oscillating mirror. The light sheet, while scanning, passed through a second cylindrical lens with the same orientation as the first to further expand the sheet. Ultimately, this approach produced slowly diverging rather than parallel light sheets (as a small angular adjustment was needed to sweep the whole volume of interest) across the volume. A similar approach was demonstrated by Lawson and Dawson (2014) but with the second cylindrical lens set at 90° to make the scanned sheets parallel. Unfortunately, approaches like this are typically limited in spatial extent as they require large, high-quality optics so volumes are typically only a few centimeters deep. The new scanning method detailed in this article allows for easy adjustment of the scanning volume, requiring no mechanical adjustment, as well as capturing data over large (0(10 −1 ) m) deep scan volumes.

Density measurements
Density fields are often measured using planar-laser-induced fluorescence (PLIF). In a standard PLIF setup, a fluorescent dye is added to passively tag one of the solutions and, through careful calibration, the concentration field of the fluoresced dye can be related to the underlying density field (Crimaldi 2008).
Simultaneous velocity and density measurements can be achieved through a careful choice of fluorescent dye, selecting the absorption wavelength of the dye to the wavelength of the laser, and using a dye with a large Stokes shift, such that the fluoresced wavelength is significantly larger than the absorption wavelength (e.g. Webster et al (2001)).
As PLIF only requires a single image for each plane, it is relatively straightforward to extend it to 3D using the scanning approach. Numerous studies have made such measurements (e.g. Dahm et al (1991) and Tian and Roberts (2003)) but, to date, simultaneous velocity and density measurements in a volume have typically been limited to small spatial extents (e.g. Krug et al (2014)). As it is desirable to have simultaneous density and velocity measurements over large volumes, we have extended the scanning approach to allow fast scanning of large volumes that readily allow the simultaneous capture of 3D density fields.

Hardware
In the present study, a dual-cavity, pulsed laser is used as the light source. The laser has a maximum repetition rate of f laser for each cavity. To accomplish the scanning, a mirror and the light-sheet-producing optics, consisting of a set of diverging (plano-concave) cylindrical lenses, are positioned on a linear traverse. The layout of all of the components that produce and scan the light sheet is shown in figure 1.
To allow overlapping light sheets, a pair of oscillating mirrors is used to superimpose a small additional translation of the light sheet in the scanning direction. To accomplish this, the beam from the pulsed laser is first directed by a fixed 45° mirror onto the first oscillating mirror. The beam is reflected by this mirror onto the second oscillating mirror before being redirected to the traverse carriage by a second fixed 45° mirror (see figure 2(a)). To facilitate initial beam alignment, the two fixed 45° mirrors are housed on adjustable mounts. Once on the traverse carriage, the beam is reflected by a final 45° mirror that passes the unexpanded beam through a set of diverging cylindrical lenses to form a light sheet perpendicular to the traversing direction (see figures 1 and 2(b)). Figure 3 shows the traverse position (black line) and light sheet position (filled circles) after taking into account the displacement induced by the oscillating mirrors. The grey line indicates the approximate (due to the finite response time of the mirrors) position of the light sheet if it were a continuous light source. The figure demonstrates how the oscillating mirrors effectively allow the continuous motion of the traverse to be broken down into a series of discrete pairs, thus allowing overlapping light sheets for the PIV analysis.
By positioning the oscillating mirrors as shown in figure 2(a), and rotating them in tandem, the beam is displaced vertically before entering the light-sheet-producing optics. The parallel displacement in the vertical, after passing through the 45° mirror at the base of the optics (see figure 2(b)), results in a parallel displacement in the traversing direction. To preserve the parallel nature of the beam displacement, the thickness of the light sheet δz should be controlled using a focusing module (e.g. a set of spherical lenses) between the laser and the oscillating mirror assembly.
For a given oscillating mirror assembly (i.e. a given L c and H c the horizontal and vertical distance between the mirror centres, respectively) the resting state angle is prescribed by θ 0 = 2θ m = tan −1 (H c /L c ). The vertical beam displacement D m /2 can then be determined from Here, θ is the angular perturbation away from the resting state θ 0 (see figure 4). The apertures of the two oscillating mirrors do not need to be equal and, as clear from figures 2(a) and 4, the aperture of the first mirror only needs to be large enough to (2) beam tube to enclose the laser beam; (3) housing containing the oscillating mirror assembly that allows the beam to be displaced vertically; (4) linear guide rail that allows the light sheet-producing-optics to travel in the scanning direction; (5) carriage containing a 45° mirror and the light-sheet-producing optics; (6) light sheet; and (7) motor used to move the traverse. (1) first fixed 45° mirror that directs the incoming beam from the laser into the plane of the sketch; (2) and (3) oscillating mirrors that allow the beam to be displaced vertically a distance D m ; (4) second fixed 45° mirror that redirects the beam down to the scanning optics. (b) Beam path after exiting the oscillating mirror assembly, looking from the side of (a): (5) beam after exiting the mirror assembly in (a); (6) fixed 45° mirror on the traverse carriage that makes the laser beam vertical (and the vertical beam displacement horizontal); (7) cylindrical lenses to produces a light sheet; (8) 3D view of (7) to illustrate the curvature of the cylindrical lens. The position of the laser beam and the carriage are shown for two sequential laser pulses, with the first and second pulses of the pair shown by the solid and dotted lines, respectively.
accommodate the incoming beam. The aperture of the second mirror A m needs to be larger (in one direction) than the beam diameter to accommodate the beam displacement, and its size (along with the resting angle θ 0 ) will determine the maximum beam displacement possible. Larger mirrors require larger galvanometers and generally have a slower response. Increasing the distance between the two mirrors allows for a smaller θ (and faster response) for a given D m , but with a larger positioning error due to the angular resolution of the galvanometer.
To simultaneously obtain the full velocity field and density field, a set of three cameras are required for the measurements (a typical layout is shown in figure 5). For the present discussion, we keep the cameras stationary as this removes a potential source of noise from the velocity fields. As is the case for non-scanned volumetric methods, having the cameras stationary requires that the lenses have a sufficiently small aperture (large f-number) so that the particles remain adequately focused across the volume of interest. For the stereo PIV measurements, two cameras are typically positioned with some angular separation and fitted with Scheimpflug adapters to provide better focusing on the scanned planes. These cameras need to be fitted with shortpass filters (or, preferably, bandpass filters centred at the wavelength of the laser) to eliminate the signal fluoresced by the dye. For the density measurements, the third camera is ideally positioned with its optical axis normal to the light sheet (to minimise distortion) and needs to be fitted with a longpass filter to remove the light directly scattered from the particles, leaving only the signal from the fluorescent dye.
The pulsed laser, oscillating mirrors, traverse and cameras are all triggered using a hardware-based timing system with all of the components sharing a common clock to keep everything synchronised (see the appendix A for more information).
To perform the PIV analysis, an accurate coordinate calibration procedure is required. This is especially true given  (2) the light sheet that is scanned a distance W; (3) the array of cameras: Cam A and Cam B for particle images and Cam C for the concentration of a scalar field for PLIF. the scanned nature of the measurements, as the optical path between the light sheet and the camera varies in time. The calibration procedure used for the measurements presented is detailed in appendix B.

Experimental setup
The methodology outlined in this paper was used to perform volumetric measurements of the buoyancy-driven exchange flow in an inclined duct similar to that studied by Meyer and Linden (2014), Lefauve et al (2018) and Lefauve (2018). We opted to use this particular experiment with our new measurement method as it is stratified (hence knowledge of the density field was desired) and has 3D spatial structure (due to the confining duct boundaries). Moreover, turbulent flow states are possible that result in a complex 3D structure of both the velocity field and density field. To this end, we chose our parameters such that we were in the 'Turbulent Regime' identified by Meyer and Linden (2014).

Experimental apparatus
The experimental setup is sketched in figure 6 and consisted of two reservoirs connected only through a duct. Each reservoir had dimensions 1000 × 200 × 500 mm and the duct had a square cross-section of height and width H = 45 mm with a total length L = 1350 mm. The duct and the reservoirs were made of Perspex (acrylic) and were of good optical quality. The duct passed through a flexible gasket that was located in the central wall separating the two reservoirs. This allowed the duct to be tilted whilst still maintaining a seal between the reservoirs. The reservoirs were filled with fluids of different densities ρ 0 ± ∆ρ/2 where ∆ρ = 11.7 kg m −3 . The duct connecting the two reservoirs was tilted at an angle θ = 5 • from the horizontal (figure 6).
We chose to orientate the x axis with the streamwise direction of the flow, along the duct, with z in the spanwise direction (the scanning direction in this setup). The y axis was then inclined at an angle of 5 • from the vertical upwards direction. Note that, rather than using the convention of having the z axis aligned with gravity as is common in stratified flow literature, our chosen orientation matches that used in the methodology of this paper (i.e. z in the scanning direction). This coordinate system had its origin in the middle of the duct, such that −L/2 x L/2 and −H/2 y, z H/2. Our measurement volume was approximately 6H (300 mm) in the streamwise direction and spanned the full vertical and spanwise extents of the duct (both H or 45 mm here). The measurement volume was offset from the centre of the duct (in the relatively dense reservoir, as shown in figure 6) to avoid the flexible gasket that separated the two reservoirs.

Optical components
For the experiments presented in this paper, a frequency-doubled dual-cavity Litron Nano L100 Nd:YAG laser was used as the light source, with a wavelength of 532 nm. The laser had a maximum repetition rate f laser = 100 Hz for each of the two laser cavities, each outputting 50 mJ per pulse. For the galvanometer mirrors we used two Thorlabs GVS311/M units (2) inclined duct that confines the exchange flow between the two reservoirs; (3) reservoir containing the relatively dense fluid; (4) reservoir containing the relatively light fluid (tagged with Rhodamine 6G for PLIF). The laser beam was emitted from the scanning system (5) and illuminated the flow through the base of the duct. Measurement volume: a 3D view of the measurement volume (1) is shown clarifying the orientation of the coordinate system used.
attached to a Thorlabs GPS011 driver/power supply. The aperture of the mirrors was A m = 10 mm and the orientations were controlled by an analogue signal that determined the subsequent rotation of the galvanometer. Given that only small rotations of the mirror (<1°) were needed in the current system, the highest resolution setting for the controller was used (a ±10 V input signal corresponded to a ±10° rotation). In the arrangement here, the mirrors were mounted so that their resting state was θ 0 = 19.9 • with a horizontal and vertical distance between the centres of L c = 188 mm and H c = 68 mm, respectively. Given this geometric arrangement, the aperture of the mirrors, and the diameter of the beam from the laser (∼4 mm), the beam could be displaced by ±2.9 mm (with a relative rotation of ±0.4° to the mirrors resting angle θ m ), which gives a maximum displacement between image pairs of D m = 5.8 mm. The angular resolution of the galvanometer introduces negligible error for our geometric arrangement, and is 0.2% of the maximal rotation ±0.4° required to have the maximal beam displacement. Note that this maximum allowable beam displacement D m (see figure 2) set an upper limit on the speed of traversing V for a fixed inter frame time ∆t as, to get overlapping light sheets, we require V < Dm ∆t . For the current system, in single-pulse mode with ∆t = 10 ms, this meant that the oscillating mirrors limited the maximum speed of the traverse to V max = 580 mm s −1 . Note that, in practice, the limit was less than this due to the limited aperture of the sheet-producing optics and the non-negligible divergence of the laser source.
To produce a scanning light sheet from the nominally circular beam emitted from the laser, a system of cylindrical lenses was mounted on the traversing carriage 1 . The traverse carriage was mounted on an Igus DryLin SAW rail system. The system was chosen due to its sliding plastic bearings, avoiding a metal on metal contact that is prone to introducing vibrations into the system, and as it was impervious to salt and water. This rail was also chosen for its relatively wide separation between the rails that minimised any roll of the carriage during motion. To move the carriage, a stepper motor was chosen for the simplicity of positional control without relying on a separate position resolver (as, for example, would be the case with a conventional servo motor). The traverse carriage was attached to a stainless steel lead screw of pitch 2 mm, diameter 10 mm, and with a maximum travel of 500 mm. Due to the inertia of the carriage and friction from the rails, the carriage could not simply be put instantaneously into motion at a constant velocity and so a constant acceleration/deceleration phase was implemented. For the experimental results shown here, the total time spent accelerating/decelerating was ∼15% of the scanning period. As a consequence, the spacing of measurements in z was non-uniform at the beginning and end of the scan.
Inevitably, vibrations were introduced into the traverse system. However, the vibrations introduced into the carriage because of the stepper motor were found to be negligible in practice (micro stepping was used to help minimise this). Moreover, to check if the traverse system deteriorated over time, a three-axis accelerometer was attached to the traverse carriage to monitor the vibration levels during the scanning. Even after multiple uses in the lab environment, the accelerometer data showed that the vibration levels of the carriage had not increased. To minimise any residual vibrational feedback from the complete traversing system onto the experiment itself, the traversing system and laser were mounted on a large (0.6 m × 1.2 m) honeycomb optical breadboard (Thorlabs PBG51507).
Three Teledyne Dalsa Falcon2 8M cameras were used to image the flow. These ten-tap CameraLink CMOS cameras had a maximum resolution of 3320 × 2502 pixels. However, given the aspect ratio of the duct, a reduced resolution for each of the cameras (3320 × 824 pixels) was used. This reduced resolution also allowed for the 200 frames per second needed to match the maximum speed of the lasers in double-pulse mode. The three cameras (two for PIV and one for PLIF) were fixed in position to one side of the duct (as shown in figure 6), and the angular offset between the two PIV cameras was chosen to be ∼80°. To improve focusing across the image, both PIV cameras were fitted with Scheimpflug adapters. It was not deemed necessary to install liquid filled prisms in between the stereo PIV cameras and the volume of interest but we note that in some situations this could be advantageous, as discussed by Prasad and Jensen (1995). The PLIF camera had its optical axis normal to the light sheet and so, unlike the PIV cameras, did not require a Scheimpflug adapter. The PIV cameras were fitted with Micro Nikkor 60 mm f/2.8D lenses at aperture f/8 and the PLIF camera had an Nikkor 50 mm f/1.2D lens at aperture f/1.2 (so that only a low concentration of fluorescent dye was required in the dyed layer, consistent with the linear relation assumed in the calibration calculation, see appendix B.4). Locating the PIV cameras 0.6 m away from the measurement volume provided an adequate field of view (the width of the section imaged was approximately 250 mm) and allowed the images to be adequately focused across the scan. Due to the increased aperture, the PLIF camera was positioned 1 m away from the measurement volume (the width of the section imaged was approximately 320 mm for the PLIF camera).
As it was not feasible to only seed particles in the region of interest, a transparent Perspex (acrylic) box filled with water (but with no particles) was positioned between the inner wall of the reservoir and the outer boundary of the duct itself in the optical path between the cameras and the light sheet. This 'optical box' removed a large portion of noise present in the raw images for all cameras that would otherwise occur due to the optical path between the light sheet and cameras encounter ing a large number of particles when high seeding densities were used.

Experimental protocol
Prior to the experiment, the system was calibrated as discussed in appendix B.1. The calibration refinement step was typically performed with particle images captured during the experiment with a second set of refinement images captured between Cam C (the PLIF camera with the filter removed) and Cam A either before or after the experiment. For the experiment, both of the reservoirs were filled with fresh water and allowed to settle overnight to help degas the fluid, to minimise bubble formation during the experiment, and reach the ambient temper ature of the laboratory. The duct was then inclined to the desired angle (as determined by a Digi-Pas DWL-280Pro digital inclinometer) and the end open to the reservoir that would contain the denser fluid was temporarily sealed (to avoid unwanted exchange of fluid between the reservoirs before the start of the experiment).
The density difference between the reservoirs was achieved using two salt solutions, NaNO 3 and NaCl, with their densities measured at a temperature of 20 °C (the same temperature as the laboratory environment) using an Anton Paar DMA 5000 density meter. This particular combination of salts was selected to allow matching of refractive index and as they, to a good approximation, mix linearly and have similar diffusivities at the low concentrations required here (Olsthoorn and Dalziel 2017). Although the heat of dilution of NaNO 3 is significantly greater than that of NaCl, for the concentrations used here the temperature change was not dynamically significant to the mixing process, and the associated refractive index change was small enough to not cause an issue. We matched the refractive indices of the solutions at 532 nm (the wavelength of the laser source) with a relative error of ∆n/n ≈ 10 −4 and verified this by a handheld refractometer (Reichert Technologies Goldberg).
For PIV, polyamide particles with a diameter of 50 µm and density between 1020-1030 kg m −3 were used to seed the flow. The particles were chosen due to their small ratio of settling to mean flow velocities V p /( √ g H) = 4.17 × 10 −4 while still being large enough to provide a clear particle image. The polyamide particles were added to the flow with a small amount (∼5 ml) of Magnum TM Rinse Aid to prevent aggregation of particles. Enough seeding particles were added to the flow to ensure that that there were always 5 particles in any interrogation window.
At this stage, the first sequence of calibration images for the PLIF measurement could be captured, as discussed in appendix B.4. Note that before recording any calibration images for the PLIF measurements, or beginning capture during the experiment, the laser and associated optical components were allowed to warm up for at least 30 s. This avoided a non-trivial time-varying change in spatial structure and magnitude of the light sheet. Such changes would otherwise hamper the PLIF calculation. After this first sequence of calibration images had been captured, rhodamine 6G was added to the reservoir containing the relatively light NaCl solution and mixed giving a final concentration of rhodamine 6G in the reservoir of C 1 = 15 µg l −1 . Adding the dye to this reservoir avoided having the light sheet pass through a layer of dyed fluid before entering the measurement volume. After ensuring the dyed fluid was mixed into the reservoir and the duct, the second sequence of PLIF calibration images were captured.
The experiment could then be started by simply removing the temporary seal from the end of the duct. Initially, there was a transient stage of the experiment as a gravity current of denser un-dyed NaNO 3 fluid propagated through the duct into the opposite reservoir. After this period of transient flow had finished we began taking the 3D measurements reported here.

Measurement resolution and processing
The in-plane resolution of the measurements for the setup discussed here were primarily set by the laser system and the choice of camera. The resolution in the scanning direction was predominantly set by the laser as the maximum number of velocity fields that could be obtained in 1 s was 100 (achieved by firing the two cavities of the laser out of phase at their maximum repetition rate 100 Hz). Moreover, the precise choice of laser and optical components used to produce the light sheet controlled the thickness of each plane δz, and the subvolume over which each of the initial 2D pixel displacements were determined. For the measurements herein, given the small depth required to image the duct (45 mm) and the typical thickness of the light sheet in the current setup (1-2 mm), we chose to discretise the domain into ∼40 planes. This avoided unnecessary overlapping of light sheets given that the uncertainty in the position of the velocity fields in the z direction will be the same order as the light sheet thickness. To adequately resolve the particle displacements, the system was in doublepulse mode so there was a spacing between light sheets of ∆t = 7 ms between the pairs of analysed frames.
The resolution of the in-plane measurements in (x, y ) was set by the resolution of the camera, the seeding density of the PIV particles, and the PIV algorithm used. For the results shown here, all of the raw images were processed using DigiFlow (Dalziel et al 2007). The processing used the DigiFlow 2017a PIV algorithm, selecting an initial interrogation window of height and width 31 pixels and a spacing of 12 pixels (both horizontally and vertically) equivalent to a 60% overlap of interrogation windows. The algorithm underlying these PIV calculations has some important differences from most PIV implementations. First, as introduced for synthetic schlieren by Dalziel et al (2000), the pattern matching kernel is based on an L 1 norm measure of the differences between images in the interrogation windows and is used in place of the more common L 2 norm of a cross-correlation function. This kernel is applied iteratively, utilising displacement information from previous passes to advect the image pair captured at t = t i ± ∆t/2 to their anticipated state at t = t i . The size, shape, and weighting profile of each interrogation window are adjusted (based on their information content and the spatial gradients in the displacement field) during this process to faithfully capture high gradients while achieving a low level of noise. Additionally, a strategy of weighting or removing anomalous pixels increases further the robustness of the results to particles entering or leaving the light sheet.
The effective resolution is not constant across the field of view due to the distortion across the images, caused by the angles of the cameras relative to the light sheet. The resolution also varies as the light sheet is scanned, a consequence of the cameras being stationary, so there is higher resolution when the light sheet is closest to the cameras. However, with the setup detailed here, the loss of resolution because of the scanning was negligible (a change in pixels of <5% for a given physical length over the depth of the scan) compared to the distortion of the images (a change of ∼20% across the image). Therefore, the maximum resolution of the velocity measurements in the (x, y) plane reported here was one velocity vector every 0.51 mm and the minimum resolution was one velocity vector every 0.63 mm. The total yield of vectors in a scan, acquired every 0.77 s, was approximately 400 × 100 × 40 (x, y, z). Indicative of the quality of the vector fields, the reconstruction error field E(x) (see equation (B.12)) had an average in the (x, y) plane of E x,y < 6.7 × 10 −2 (spanning all z in the scan) with a maximum standard deviation σ max = 4.6 × 10 −2 , both in units of pixel/frame. It was observed that the error systemically increased in z (increasing distance from the cameras). However, it is worth noting that, due to size constraints of the current dual-plane calibration target, the calibration could only be performed over half of the spanwise extent of the duct (the side situated closest to the cameras). This is a possible cause of the increasing trend observed as a similar trend is also found if we simply go between the forward and inverse mappings as defined in appendix B.2.
The density fields in world coordinates were calculated as outlined in appendix B.4. Note that, for the data shown here, only the first of the two PLIF images was used to reconstruct the density field. As no interrogation window is required for the PLIF calculation, the resolution of the PLIF measurements is higher than those of the velocities. In the current setup, there was a density measurement in the (x, y) plane every 0.1 mm. In total, there were approximately 3000 × 500 × 40 (x, y, z) density measurements every 0.77 s, given by the resolution of the camera and the number of frames in a scan. For both the velocity and density fields, the resolution in the scanning direction z is 1.26 mm with an accuracy set by the light sheet thickness δz (1-2 mm).
Before presenting the data, the density fields were processed to remove line artifacts from (x, y) planes (due to light rays passing through residual air bubbles, clusters of particles, and imperfections in the tank surfaces) using a threestep approach. As all of the rays emanated from some virtual origin below the tank and spread with the light sheet, the first step straightened the lines by mapping the (x, y) planes into a 'ray coordinate' system producing planes with vertical line artifacts. The second step removed these vertical line artifacts from the planes using a wavelet method described by Münch et al (2009). The planes, now without line artifacts, were mapped back to world coordinates. Finally, after removing the line artifacts, the density data was interpolated onto the lower resolution grid of the velocity data after being median filtered over a suitably sized window.

Experimental results
For the results presented here, we choose to non-dimensionalise velocities u by √ g H, where g = g∆ρ/ρ 0 is the reduced gravity, and normalise all lengths by H/2. As a result, with the duct angle fixed at θ = 5 • , we can construct a Reynolds number Re = √ g HH/(2ν) = 1516 to characterise the flow. The natural time scale to non-dimensionalise time is the advective time 1 2 H/g . Finally, the non-dimensional density field is ρ = 2(ρ − ρ 0 )/∆ρ such that −1 ρ 1. Therefore, our final measurement volume in non-dimensional units was −17.4 x −6.3, −1 ỹ 1, −1 z 1 and a total of 561 advective time units were captured with a non-dimensional time between volumes of 2.40. The tilde indicating non-dimensional variables will be dropped henceforth.
To get an overall impression of the flow, time-averaged quantities of u and ρ are shown in figure 7. In comparison to a laminar exchange flow, where there are two layers separated by a sharp density interface, mixing is evident in this relatively high Re number flow. Here, thin and approximately well-mixed layers are confined to the horizontal boundaries at y = ±1. These layers are separated by a weakly-stratified interior (figures 7(a) and (d)). Mean profiles of the streamwise velocity u are also shown in figures 7(b) and (h) and indicate the spanwise variation present within the flow due to the confining lateral boundaries at z = ±1. In the same geometry, the importance of this confinement to the flow at lower Re has been investigated by Lefauve et al (2018), who found a non-trivial modification to the classical Holmboe instability with significant 3D structure. These results demonstrate that, even at relatively low Re, 3D spatial measurements provide meaningful physical insight that would be lost in single plane measurements. Figure 7(g) also shows that there is no significant spanwise variation of the density, and therefore no mean pressure gradients in the spanwise direction due to horizontal density gradients. Finally, indicative of the stability of the flow, the mean gradient Richardson number Ri g (y, z) is also shown in figures 7(c), (f), and (i). Note that here where is the buoyancy frequency associated with the averaged density profile ρ x,t and is the squared shear associated with the mean profiles of streamwise velocity u x,t , with · x,t indicating the average in the streamwise direction x and time. From figures 7(c) and (f) it is evident that there is an approximately constant region where Ri g (y, z) 0.12 that is confined vertically to the centre of the duct, associated with the weakly-stratified interior and approximately constant shear in this region (see figures 7(d) and (e)). A low value of Ri g (y, z) is expected given the turbulent nature of the flow, where values are typically Ri g (y, z) < 0.25 (Holt et al 1992). Instantaneous data are shown in figures 8 and 9. Significant 3D variation is evident in all components of velocity u throughout the entire domain, as shown in figures 8(a)-(h). At the same instant in time and on the same planes, the density field ρ and the enstrophy |ω| 2 are shown in figures 9(a)-(e). These figures highlight the necessity of time-resolved 3D fields in analysing such flows, especially if coherent structures are to be extracted from the complex flowfield. Furthermore, the enstrophy could not be calculated without knowing gradient information of all components of u.
To demonstrate the near-instantaneous nature of the volumetric data, isosurfaces of enstrophy and density are shown in figure 10. To be able to capture structures, we only need to resolve the Eulerian timescales of the flow. This is distinct from needing to resolve the Lagrangian timescale between pairs of images to accurately determine the particles displacements (which are Lagrangian tracers) that would require a much faster frequency of acquisition over each volume (for example, as required in conventional non-scanned volumetric methods). Crucially, quantities in the scan direction (including derivatives) are reliably resolved as the time between velocity/ density subvolumes are only separated by 2∆t (therefore only double the ∆t needed to resolve the particle displacements) and the data are effectively skewed across the scan direction. The distorted nature of the data could be advected using the velocity information but this step has not been carried out on the data presented here.
A further measure of the quality of the data can be inferred from ζ = ||u − u div || 1 the L 1 norm of the difference between u and a divergence free field u div calculated using the method described by Wang et al (2017). For the data shown here, the mean of this positive quantity over the whole volume and all time is ζ = 2.5 × 10 −2 with standard deviation σ = 2.2 × 10 −2 (all in non-dimensional units) illustrating that the correction is small (note that the value of the mean corresponds to ∼0.17 pixel/frame, the same order as our error estimate in section 4.4). Given that all experimental data will have some non-zero divergence, and the fact that the measurements are near-instantaneous rather than truly instantaneous, we anticipate that ζ > 0. However, our measurements are validated as only a small numerical correction is required and we note that the qualitative difference between the u fields shown in this paper and the corresponding u div fields is almost imperceptible.

Discussion and conclusions
Depending on the flow of interest, there are two enhancements that extend the methodology presented in this paper. These enhancements make use of the same fundamental approach but facilitate its use with 'slower flows' (flows that are slow compared to the camera frame rate) and flows with a strong out-of-plane motion.
For flows that are slow relative to the camera frame rate, the methodology outlined so far would be slaved to the time between light sheets. Slow flows require a larger ∆t so that the particle displacements are sufficient to obtain accurate velocity measurements. Therefore, for a given number of light sheets in a scan N the time taken to scan the volume will be N∆t. However, so far we have only discussed the simple 'mode 0' operation of the system, where the mode number is given by m in z k = Z 2k = Z 2(k+m)+1 .
(6) Here, z k defines the z position where two light sheets overlap, with Z k the z position of the kth laser pulse. Therefore, for mode 0, overlapping pairs of images each comprise of an even-numbered frame and its immediately following odd number, e.g. z 0 = Z 0 = Z 1 , z 1 = Z 2 = Z 3 , etc. However, for slower flows it is beneficial to operate in higher modes where the images that are spatially coincident are separated in frame number. Operating in higher modes allows for greater temporal resolution for the complete scan that would otherwise be slave to the ∆t required to obtain a sufficient particle displacement between images. Essentially, the effective ∆t used in the PIV calculation can be increased without having to decrease the scan rate or the number of subvolumes the volume is discretised into. However, operating in higher modes does make it necessary to displace the laser beam further and therefore requires bigger oscillating mirrors and larger aperture cylindrical lenses to accommodate the larger amplitude beam displacement.
In traditional PIV methods (planar or stereo), strong outof-plane motion can cause errors in the PIV algorithm (due to Figure 10. Instantaneous isosurfaces of (a) ρ = 0.5 and (b) |ω| 2 = 12 in (b). For the data shown here, an isotropic 3D Gaussian filter with σ filt = 1 vector spacing was used to smooth the data with a central difference scheme used to compute the derivatives of the velocity field. For both plots, the vertical y axis has been stretched by a factor of two to aid visualisation and the vertical and spanwise extents of the data have been limited to [−0.75, 0.75] to avoid the signal at the boundaries obscuring the view. substantial loss of particles between frames) or demand a high acquisition rate (to minimise the loss of particles between frames) that increases the noise of the measurements due to the inevitable lower in-plane pixel displacement. A potential use of the system, be it scanning or not, is the use of the mirrors on galvanometers to allow thin light sheets in configurations that would otherwise require a thick light sheet to accommodate the out-of-plane motion w. If it known beforehand that the flow to be measured has a unidirectional mean flow in the out-of-plane direction w, the mirrors can be used to displace the light sheet in the z direction a suitable amount to accommodate the motion, i.e. D m ∼ w∆t. The anticipated out-of-plane motion could be determined by an initial experiment, with a thick light sheet, and then repeated with a thinner light sheet to improve the accuracy of the measurements.
To conclude, we have presented a highly versatile scanning PIV/PLIF technique. Our approach allows for fast, repeated scans of a measurement volume. The volume is assembled from a series of subvolumes that allow for the simultaneous measurement of all three-components of the velocity field and the density field. The size of the plane (x, y ) can be readily adjusted by simply changing the cylindrical lenses, and the distance in the scanning direction can be adjusted 'on-the-fly' by suitable commands to the traverse. Moreover, the system lends itself to complex scan sequences, e.g. varying volume sizes in the scanning direction, or several scans followed by fixed plane measurements with higher temporal resolution.
The novel addition to the setup is the use of two mirrors on galvanometers to position the light sheet during the scanning. These mirrors allow thin, overlapping light sheets between pairs of images despite the light-sheet-producing optics being continuously translated by a traverse system to scan a large volume.
Ultimately, the setup detailed in this paper can scan a volume of approximately 0.1 m 3 in less than a second, where the time taken is limited by the repetition rate of the laser and the desired number of frames in the scanning direction. Within each scanned volume, the total yield of velocity vectors is O(10 6 ) with O(10 7 ) simultaneous density measurements.
We have demonstrated the capability of our technique by taking measurements of a stratified shear flow. As the purpose of this article was to describe the method, we refer the reader to Lefauve et al (2018) and Lefauve (2018) for further details of the dynamics of the flow. Measurements of this nature are useful in the investigations of such flows because of their inherent 3D nature, due to turbulence, and the importance of the density field in such problems.
target (a plane of dots or patterns with known spacing) to several z positions, or by using a multi-plane target where the target has two (or more) planes of dots or patterns at known z locations (Prasad 2000). We opt to use a dual-plane target herein, and so, by imaging the calibration target within the experimental setup, 2D mapping functions G k± i can be calculated. The mappings G k± i map from the pixel coordinates of each camera (B.2) Figure A1. Schematic of the timing sequence sent to the various components of the system when in single-pulse mode. The laser pulse is triggered from the rising edge (+ve) of the pulse train to the lasers. An image pair (one plane of velocity) is indicated by the shaded region. Figure A3. Schematic showing the timing signal over the course of a complete scan. Acceleration/deceleration phases are shown for the traverse, along with the direction signal that controls the direction of travel for the traverse. A sync pulse is shown that is used to trigger the recording of a sequence of images for all cameras at the start a scan. Figure A2. Schematic of the timing sequence sent to the various components of the system as in figure A1 but now in double-pulse mode.
Here, P 2 ⊂ R 2 is the pixel coordinate space, on the two z planes of the target z k ± ∆z 2 , where z k is the mid-plane of the target and ∆z is the spacing between the two planes, and where W 3 ⊂ R 3 is the world coordinate space. Specifically, where i is the camera identifier and k± distinguishes between the two planes of the calibration target positioned at z = z k ± ∆z 2 . Thus there are two mappings, G k+ i and G k− i , for each camera. In general, the mapping functions G k± i are not known a priori and are commonly determined by least squares fitting of a polynomial function between the known world coordinates x of the features on the calibration target and the corresponding pixel coordinates of each camera X i .
For a non-scanning system, with a light sheet centred at z k , the geometric mapping G k i at z k can be inferred by linear interpolation between the two mappings at z k ± ∆z 2 : To calculate the world displacements ∆x, we need a 3D mapping. We opt to use an approach similar to Soloff et al (1997) and start by imaging the dual-plane calibration target with both PIV cameras. Then, by identifying a common point on the calibration target, the four pixel coordinates (two from each camera) are used to determine the 3D mapping where the subscripts A and B distinguish between the two PIV cameras, and k defines the mid-plane location of the calibration target z k where the mapping is calculated. Given the two planes (separated by ∆z) of the dual-plane target, a polynomial function can be fitted, for a given target position z k , with high order dependence in x and y (dependent on the number of features of the target), and linear dependence in z. Given the thin light sheet, higher order z dependence is generally not required, but could be achieved by including more planes on the multi-plane target or by accurately traversing a planar target in z. The world velocities u = (u, v, w) on a plane positioned at z = z k can then be found from where U i (X i , z k ) = ∆Xi ∆t are the 2D (pixel) velocities for camera i and J k i = ∂F k ∂Xi is the (3 × 2) Jacobian matrix (with units of world/pixel) associated with the mappings at the z k location of camera i. For the scanning system, the mappings and Jacobian matrices are required for all z in the scan. To this end, the dual-plane target is positioned at a number of z k locations within the volume of interest. This enables the mappings G k i and Jacobians J k i (plus inverses of both) to be calculated at several z locations spanning the volume. A least squares fit in z of the polynomial coefficients of the mapping functions is then used to generate z-dependent 2D mappings G i (X i , z) and z-dependent Jacobians J i (X i , z) over the whole volume to be scanned.
The first step in calculating stereo velocities is calculating the 2D velocities in pixel coordinate space P 2 for each of the PIV cameras at a given z location U i (X i , z). The velocities U i (X i , z) in pixel coordinate space P 2 at a known z location are mapped to a common grid in world coordinate space W 3 . The z position of the velocity fields is used so that the pixel velocities of the two cameras in pixel coordinate space P 2 can be mapped to the corresponding world coordinate space W 3 using the z dependent 2D mappings generated during the calibration: The Jacobian J i , initially in P 2 space separately for each camera, is also mapped to world coordinate space W 3 in the same manner: We do this for each pair of 2D2C planes produced by the 2D PIV algorithm for cameras A and B to construct a sequence of 2D3C planes at different z locations in the scan. These 2D3C planes are finally combined to construct volumetric 3D3C measurements.

B.2. Error
For each z location of the scan, an estimate of the error in the stereo reconstruction can be determined by back-projecting the world velocities û onto the velocities (in pixels/frame) of the two cameras. For convenience, we do this calculation on the common world grid to obtain Û * i as follows: where Ĵ −1 (x) is the inverse Jacobian in world coordinate space W 3 found from is the (4 × 3) Jacobian matrix (with units of pixel/world) of the inverse 3D mapping F −1 i . A measure of the error can be obtained by comparing the back-projected velocities with the velocities calculated by the PIV algorithm mapped to W 3 . Assuming equal weighting among each of the velocity components, we construct the error field from the magnitude of the vector difference: This error estimate can be used as an additional quality check of the velocity fields by removing vectors where the reconstruction error in E(x) exceeds some threshold value (typically chosen to be 0.5 pixels/frame as in Wieneke (2005)). An estimate of the error introduced by the mappings can also be found using the forward mapping on the back-projected fields Û * i to obtain û * i and then the inverse mappings on this field to obtain a second set of back-projected mappings. These back-projected mappings can then be compared to Û * i to highlight possible errors associated with the mappings.

B.3. Calibration refinement
An important step in calculating the velocities using stereo PIV is the calibration refinement technique. The stereo reconstruction of velocity, as outlined in appendix B.1, relies on the calibration between each camera and the world coordinates being accurate. This ensures that a velocity perceived from one camera is correctly reconstructed with a velocity from the other camera, i.e. the two velocities are coincident in world coordinates. Although, in practice, the calibration procedure is performed carefully, it is difficult to perfectly align the calibration target with the position of the light sheet and even small discrepancies can yield large errors in the stereo reconstruction (Wieneke 2005). With the scanning system, this is complicated further as any small misalignment between the laser beam and the optics on the traverse carriage can lead to a systematic positional error in the light sheet. To correct for this, we opt to use an iterative two-step calibration refinement technique using a methodology similar to that of Willert (1997) for a fixed light sheet.
Effectively, an algorithm similar to PIV is run on synchronised images from the different cameras all projected on to world coordinates. In the first step, the mean disparity between the images is used to refine the light sheet position, i.e. the position of each pair of light sheets is updated until the mean disparity between the images is minimised. The second step performs finer adjustment across the whole image, to effectively remove any residual skewing of the light sheet relative to the coordinate system. This ensures that the 2D velocity information from each of the PIV cameras overlap. In general, this two-step approach should be iterated but typically we find that any residual error in position is much smaller than δz after one pass.

B.4. PLIF calculation
To determine the density field, the signal from the fluorescent dye added to the flow has to be calibrated. Assuming negligible attenuation and a linear camera response (see Crimaldi (2008) for more details), the calibration of the experimental images is simply C(X, z) = I(X, z) − I 0 (X, z) I 1 (X, z) − I 0 (X, z) C 1 , (B.13) where C(X, z) is the calculated concentration field, I 1 (X, z) is a reference image containing a known homogeneous mixture of dye that can be used to remove any spatial variation present in the light sheet and establish the concentration of dye in the images, and C 1 is the dye concentration for which I(X, z) = I 1 (X, z). These two reference images are determined as follows. First, a scan-position-dependent background reference image I 0 (X, z) is calculated at each z position of the scan by recording O(10) scans of the volume containing no fluorescent dye and taking the mean at each z l location. To minimise errors, no fitting is used for the background images and a distinct image is found for every z location in the forward and backward scan sequences. In a similar manner, a second reference image I 1 (X, z) is calculated by averaging over O(10) scans of the volume containing a homogeneous mixture of the highest concentration of dye in the experiment C 1 , again producing 2N images. Finally, the concentration measurements found from equation (B.13) are mapped to world coordinate space W 3 . As discussed in appendix B.1, this is achieved in the scanning system using a least squares mapping G C (X C , z) calculated from the images recorded of the calibration target for the PLIF camera (Cam C). The concentration Ĉ in world coordinate space W 3 can then be found from (B.14) As for the PIV cameras (Cam A and B), the coordinate system of Cam C can be refined using the strategy discussed in appendix B.3. In this case, only the second step of refinement is needed, and the coordinate system of Cam C can be refined using simultaneously acquired particle images from Cam C, with its longpass filter removed, and one (or both) of the PIV cameras. To ensure all of the coordinates systems lie in the same plane, the disparity map is calculated using the refined coordinate system of Cam A or B (or both) with the full disparity map (or the average disparity map if the refinement is conducted with both PIV cameras) applied to Cam C images only.
To relate the concentration to the density field ρ we have ρ = ∆ρ C 1Ĉ + ρ min , (B.15) where ∆ρ = ρ max − ρ min and ρ max , ρ min are the densities of the fluid with the maximum concentration of dye and the fluid containing no dye, respectively. It is worth noting that, in the scanning system, we have twice as many ρ fields as velocity fields because the PLIF calculation only requires a single image (compared to the two needed to calculate a velocity field using PIV). That being said, there is a choice in how to make use of the extra information. We could use either a single ρ field, corresponding to the first or second raw particle image used in the PIV calculation, or the mean of the two ρ fields could be taken to help remove random noise from the data. Alternatively, the velocity information calculated for the corresponding pair of ρ fields could be used to advect the fields forward and backward in time by ± 1 2 u∆t before taking the mean of the two fields.