Reynolds stress tensor measurements using magnetic resonance velocimetry: expansion of the dynamic measurement range and analysis of systematic measurement errors

This study presents magnetic resonance velocimetry (MRV) Reynolds Stress measurements in a periodic hill channel with a hill Reynolds number of Re = 29,500. The velocity encoding scheme is based on the ICOSA6 method with six icosahedral encoding directions and multiple encoding values are measured to increase the dynamic range. The full Reynolds stress tensor is obtained from a voxel-wise three-dimensional Gaussian fit using the magnitude data of all acquisitions. The MRV results are compared to a wall-resolved large eddy simulation and laser Doppler velocimetry measurements conducted in the same channel. It is shown that the MRV Reynolds stress data have excellent precision and agree qualitatively with the reference data. However, there are apparent systematic deviations. One of the most prominent error contributions is the signal attenuation caused by higher orders of motion, which leads to an overestimation of the turbulence level. Another fundamental error is identified in the assumption that the turbulence is Gaussian distributed. With the presented reconstruction technique, the MRV data are fitted to a statistical model, and depending on the examined flow setup, the Gaussian model can lead to considerable errors. Possible ways of how to reduce all identified errors are presented. In summary, this technique enables Reynolds stress tensor measurements in complex internal flows with high dynamic range and excellent precision. However, several issues need to be resolved to make the turbulence quantification more accurate.


Introduction
The measurement of the Reynolds stress tensor (RST) is an important topic for the investigation of turbulent wallbounded flows, and this can be achieved with magnetic resonance velocimetry (MRV) (Haraldsson et al. 2018). MRV does not require physical or optical access to the flow field, which enables simple experiments in highly complex channel geometries. However, previous studies reported that MRV-based turbulence measurements have possible bias, which limits the use of such data (Elkins et al. 2009). The aim of this study is to analyze the sources of error in the Reynolds Stress data obtained with MRV and to further develop this measurement technique.
Accurate turbulence measurements require a measurement system that has a high dynamic velocity range and has low measurement uncertainty. Hot wire anemometry (HWA) and laser Doppler velocimetry (LDV) are both point-wise (1D) techniques and often serve as the ground truth in turbulence quantification (Shirai et al. 2006;Samie et al. 2018). HWA is well suited for velocity measurements with high dynamic range, but being an intrusive technique, it comes with several experimental challenges. LDV is a non-intrusive measurement technique but requires seeding particles in the flow. Furthermore, Reynolds Stress measurements with LDV and HWA typically result in low data rates because of coincidence filtering, i.e., the same turbulent eddy must be detected by multiple channels to obtain the Reynolds shear stresses.
Particle image velocimetry (PIV) and particle tracking velocimetry (PTV) can provide planar (2D) or volumetric (3D) information of turbulence statistics. Compared to HWA and LDV, these techniques are known to have relatively high random errors which can bias the calculation of the Reynolds Stresses (Wilson and Smith 2013). Several sources of error must be carefully compensated for in order to ensure sufficient accuracy (Scharnowski et al. 2019). In the case of internal flows, velocity measurements with PIV and PTV typically require that major parts of the channel walls are transparent, and that the refractive index of fluid and channel material is matched. Full three-dimensional Reynolds Stress measurements, which are here termed 3-dimensional 6-component (3D6C) RST measurements, are typically out of reach in these experiments.
The experimental challenges with MRV are entirely different. The working fluid must contain a measurable nuclear spin and magnetic susceptibility of the fluid and all surrounding material must not change substantially to avoid magnetic field distortions. Also, temperature can affect the magnetic properties. Most commonly, MRV experiments are comprised of an isothermal water flow through a channel made of non-metallic material. The measurable flow velocities are typically in the range of 0.01-10 m/s but higher or lower flow velocities are possible (Elkins and Alley 2007). The main advantages of MRV are that it does not require optical access or seeding particles.
When it comes to quantifying turbulence, the MRV technique can benefit from several other factors. A coincidence filtering of instantaneous velocity samples as in LDV or HWA is not necessary because the MRV data represent the velocity spectrum of all water molecules in each voxel, not individual tracers. However, there are several sources of error that are specific to MRV, and many of these contributions might not be known yet. According to the authors' knowledge, to date, no study has yet been dedicatedly conducted to investigate the measurement errors in MRV Reynolds Stress data systematically.
This study presents MRV experiments in a periodic hill channel with a hill Reynolds number of Re = 29,500. A novel routine for data acquisition and data processing is used to reduce the measurement uncertainty and to increase the dynamic range of the Reynolds Stress measurements. The MRV results are compared with a wall-resolved Large Eddy Simulation (LES) that was set up to match the experiment. Furthermore, LDV turbulence measurements serve as the ground truth. It is shown that the MRV Reynolds Stress data have excellent precision and agree qualitatively with the reference data sets. However, there are apparent systematic deviations, which can be either attributed to the MRV encoding process or they are caused by inaccurate assumptions. Alternative measurement and post-processing methods to reduce these errors are discussed.

Fundamentals and previous studies on MRV turbulence quantification
In phase-contrast-based MRV, the voxel-wise mean velocities are encoded in the signal phase via magnetic field gradients. The design parameter is the first gradient moment (m 1 ) of these gradients, which enforces a linear relationship between the flow velocity and the signal phase. With the use of Fourier velocity encoding, the velocity spectrum within a voxel can be sampled (Moran 1982;Callaghan et al. 1988;Li et al. 1994;Newling et al. 2004). However, due to the high dimensionality of the sampled space (six velocity dimensions comprised of three variance components and three covariance components) this method requires long acquisition times and may not be feasible for many applied experiments. MRV turbulence measurements are usually based on a priori assumptions on the shape of the intra-voxel velocity spectrum in order to significantly reduce the acquisition time compared to Fourier velocity encoding. Assuming a Gaussian velocity distribution, the six-dimensional space sampled by Fourier velocity encoding can be considerably reduced to six individual three-dimensional measurements that differ in the applied first gradient moment and one reference measurement (Dyverfeldt et al. 2009;Haraldsson et al. 2018). Given the analytic solution of the Fourier transform of a Gaussian, the velocity distribution can be directly obtained from the magnitude data of these seven measurements.
As the underlying principle, turbulence within a voxel leads to a range of signal phases that add up incoherently and lead to a signal attenuation, hence decreased magnitude values. Similar to velocity-encoding, the signal attenuation depends on the direction and magnitude of m 1 , which enables a quantification of the full RST. Note, that turbulence measurements in MRV share the same measurement principles as diffusion-weighted imaging and diffusion tensor imaging, as more medical-focused MR techniques (Le Bihan et al. 2001).
The dependencies between signal attenuation, sequence parameters and flow turbulence have been applied in early studies to derive a measure of the turbulence level in pipe flows and free jets (Kuethe 1989;Gao and Gore 1991). A more quantitative study was conducted by Gatenby and Gore (1994) who found that the signal attenuation also depends on the temporal autocorrelation function of the velocity changes. They introduce a weighting function that includes the Lagrangian integral timescale, which is a measure of the average duration of each velocity change. However, such timescale is typically not derived easily and requires additional measurements (Kuethe and Gao 1995). Other researchers simply overcome this problem by considering only cases in which the timescale of the fluctuations is either sufficiently small (Newling et al. 2004), or much larger with respect to the encoding time (Elkins et al. 2009;Dyerfeldt et al. 2006). Elkins et al. (2009) conclude that most experiments conducted with MRV count to the case that the Lagrangian integral timescale is much larger than the encoding time. MRV experiments are typically carried out in large-scale models in which the flow velocities rarely exceed 1 m/s. With encoding times in the order of milliseconds, signal attenuation in such flows becomes largely independent of the Lagrangian integral timescale, which simplifies the experiment tremendously.

Theory
The central assumption in the applied MRV technique is a Gaussian velocity distribution within each voxel. By further assuming homogenous flow conditions within a voxel and a sufficiently large Lagrangian integral timescale, the corresponding image magnitude S in a voxel can be expressed as: where ⃗ k v = m 1,x m 1,y m 1,z T is the three-component encoding vector, γ is the gyromagnetic ratio, Σ the variance of the Fourier transform of the three-dimensional Gaussian velocity distribution and S 0 the signal of a velocity-compensated measurement. The parameter Σ −1 yields the RST, which contains the variance (Reynolds normal stresses) and covariance (Reynolds shear stresses) of the three-component velocity fluctuations.
In this study, the RST is obtained from six velocityencoded measurements and one velocity-compensated measurement according to the ICOSA6 encoding scheme (Zwart and Pipe, 2013;Haraldsson et al. 2018): where = (1 + sqrt(5))∕2 and m enc 1 is the applied encoding moment. The velocity encoding axes correspond to the vertices of an icosahedron with radius m enc 1 . Compared to orthogonal encoding schemes, the measurement efficiency is increased since the non-orthogonal encoding directions in the ICOS6 scheme share information, which decreases noise. However, with a fixed m enc 1 , the relationship between signal attenuation and turbulence is limited by an upper and lower limit. In low-turbulent flow regions, the signal attenuation is often not measurable. Whereas, in the high-turbulent areas, the attenuated signal falls below the noise floor. The same problem exists for the quantification of the full RST. Turbulence in wall-bounded flows is often anisotropic, and not all Reynolds Stress components have a similar magnitude. Therefore, non-uniform flows with anisotropic turbulence require a high dynamic range to quantify the full tensor precisely.
Measurements with multiple m enc 1 have been conducted by others to increase the dynamic range of component-wise Reynolds Stress measurements (Elkins et al. 2009). In that approach, a one-dimensional Gaussian function is fitted onto the magnitude data obtained with different m enc 1 to obtain the Reynolds Stress component in one encoding direction. A similar modification can be realized for the ICOSA6 encoding scheme.
In this study, the dynamic range of the Reynolds Stress measurement is increased by measuring multiple icosahedral shells with different radius m enc 1 , which are chosen to cover the entire range of all Reynolds Stresses in the flow domain. The full RST is then obtained by fitting a three-dimensional Gaussian function to all m enc 1 shells by the least-squares approach. The main difference to the scheme in Elkins et al. (2009) is that the RST is obtained from a single three-dimensional data fit compared to six individual data fits, which should increase precsision as more data points contribute to the data fit.

Methods
The following paragraphs present the flow experiment (Sect. 3.1) and describe the MRV routine used to obtain 2D6C RST data (Sect. 3.2). As reference, measurements with LDV (Sect. 3.3) and a wall-resolved LES (Sect. 3.4) are conducted in the same flow geometry under identical boundary conditions.

Experimental setup
Flow measurements were performed in a periodic hill channel shown in Fig. 1. This test case was chosen because of the high non-uniformity and anisotropy of the RST tensor. The geometry of the hills is similar to the Ercoftac case 81 "Flow over periodic hills", which is described in Temmerman and Leschziner (2001) and Jang et al. (2002).
The original test geometry is a 2D problem that can be easily modeled in CFD by applying periodic boundary conditions in the span-wise direction. In the experiment, however, it would require an extremely wide channel to exclude 3D effects. For example, Rapp and Manhart (2011) chose a channel width equivalent to 6 channel heights. Such a flat rectangular channel would not fit inside the circular bore of the MRI scanner, without decreasing the channel height and losing resolution.
For this reason, it was decided to update the geometry and create a new test case. In this study, the 2D hills were placed in a square cross-section channel to intentionally incorporate 3D flow effects. The channel height and width were 74 mm. All other geometric parameters were scaled according to the original test geometry, resulting in a hill height of 24.4 mm and a span between consecutive hills of 219.4 mm.
Since the nature of the flow problem had already changed, i.e., from a 2D problem to a 3D problem, it was also decided to increase the Reynolds number from Re = 10,595 to Re 29,500. This number is based on the hill height and the axial bulk flow velocity measured at the hill crests.
The periodic hill channel consists of 15 identical consecutive hills to provide a sufficient periodicity of the flow. All measurements were performed between the 12th and 13th hill. For comparison, Almeida et al. (1993) and Rapp and Manhart (2011) used a channel with 10 consecutive hills with measurements conducted between the 6th and 8th hill. Since the true periodicity in the experiment by Almeida et al. (1993) was found questionable (Temmerman and Leschziner 2001), the total number of hills was increased here.
The channel is made of hydraulically smooth Polymethylmethacrylat sheets. The contour of the periodic hills is realized with laser-sintered parts made of polyamide powder, which were smoothed with sandpaper until the arithmetic mean surface roughness was less than 5 μm. Upstream of the channel assembly is a cylindrical settling chamber with a diameter of 200 mm and a nozzle that smoothly contracts the flow to the square cross-section of the periodic hill channel. A diffusor is placed downstream of the channel to provide neutral pressure conditions at the outlet. Furthermore, screens are placed at the inlet and outlet of the channel section to further homogenize the flow boundary conditions. The flow supply system is comprised of two 5.5 kW centrifugal pumps which run in parallel and are driven by frequency converters. A tank with 1000 l water acts as a buffer to keep the fluid temperature stable. The flow rate is monitored with an ultrasound flow rate sensor (Deltawave C-F, Systec Controls, Puchheim, Germany) with a total tolerance of 1.8 l/min or 1.5% of the measured flow rate, depending on which is higher. Also, standard pt100 temperature sensors and pressure transducers are installed in the flow loop to monitor the flow conditions. The same setup is used for the MRV and LDV experiments, including all described components. The working fluid is purified water with a small concentration of 1 g/l Copper sulfate. This contrast agent amplifies the nuclear magnetic resonance signal from the water, but at these low concentrations has no measurable effects on the fluid mechanical properties relevant to this study. After performing the MRV measurements, a small amount of Vestosint (Evonik Industries AG, Essen, Germany) with 5 μm particle size was added to the fluid as a tracer for the LDV experiments.
Before each measurement, the flow system was kept running for several hours until stable flow conditions were reached. The final flow rate and fluid temperature are 255 l/min and 22 °C, which yields the target Reynolds number of Re = 29,500. All relevant flow conditions are listed in Table 1.

MRV measurements
A 2D slice in the symmetry plane of the channel was acquired on a 3 T MRI system (Magnetom Tim Trio, Siemens, Germany) with 40 mT/m maximum gradient amplitude and 200 T/m/s maximum gradient slew rate. The field of view (FOV) of all MRV acquisitions covered the flow between the 12th and 13th hill including several centimeters up and downstream of this section. Two receive-only coils were placed below and above the measured section. The FOV and image resolution are the same for all acquisitions.
The full 6-component RST tensor was measured with a custom phase-contrast Gradient Recalled Echo (GRE) sequence, employing the gradient design proposed in Nishimura et al. (1991). More details on the sequence are provided in Schmidt et al. (2020). The velocity encoding values on all encoding axes were calculated according to the ICOSA6 encoding scheme in Eq. (2). A total of 12 m enc 1 values were measured to achieve a sufficient signal attenuation in all flow regions. The data was averaged over 256 acquisitions to reduce measurement uncertainty to a minimum. This was mainly necessary because the turbulence leads to inconsistent phase contributions for different k-space lines, which in turn manifests itself as highly pronounced ghosting artifacts in the image. Additional measurements without flow but otherwise identical settings were performed before, after and in-between the flow measurements to measure the background phase and possible drift. The main acquisition parameters is provided in Table 2.
In this study, a voxel size of 1 mm × 1 mm × 6 mm (in x, y, z) was chosen as a tradeoff between measurement time and spatial accuracy. Note that the voxel size is uniform. Unavoidably, there are voxels at the border between flow domain and wall material that contain a mixed signal of both regions. These voxels still represent the true fluid velocity if the signal from the wall material is zero, which is here the case. As a result, MRV provides velocity data in the vicinity of the walls.

MRV data processing
The image magnitude and phase were reconstructed from the multi-channel MRI raw data using coil sensitivity maps obtained from the ESPIRiT approach described in Uecker et al. (2014). For the reconstruction of the mean velocities, the phase data were corrected for background phases by the measurements without flow. The velocity vector field was then reconstructed from the data corresponding to the two lowest m enc 1 values to minimize phase bias from higher orders of motion (e.g., acceleration).
The RST was quantified by two different methods. First, by fitting a three-dimensional Gaussian distribution to the magnitude data of the ICOSA6 encoded measurements with all 12 m enc 1 values, normalized by the velocity-compensated measurements. Data points with less than 10% normalized magnitude were excluded from the fit, to avoid voxels that are too close or within the noise floor. These voxels contain rectified noise, which would lead to an underestimation of the signal attenuation. For demonstration, the RST was also obtained by calculating the three-dimensional Gaussian   2.5, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 60]  individually for every m enc 1 value, which resembles the approach in Haraldsson et al. (2018). All data reconstruction and post-processing were conducted in Matlab R2018a (The Mathworks, Natick, USA).
An important parameter to assess the convergence of the turbulence statistics is the total sampling time. The raw image data obtained with MRV was reconstructed from 120 samples (phase encoding steps), each obtained over a readout period of 2.2 ms (reciprocal value of the readout bandwidth provided in Table 2). Considering that 256 measurements were averaged, the sampling time for each encoding point yields 67.6 s. The RST data includes the data from six encoding directions. Depending on the position in the flow field, a minimum of 33 MRI raw data sets was used for the calculation of the RST data, resulting in a minimum sampling time of 2230 s. For comparison, the characteristic convective timescale in the periodic flow is 0.28 s. This number describes the ratio of the distance between two consecutive hills and the bulk flow velocity in the unobstructed channel. As a result, the total sampling time of the MRV measurement is equivalent to a minimum of 7964 convective times, which can be considered more than adequate to achieve a sufficient convergence of the turbulence statistics.

LDV measurements
Using the same flow system as in the MRV experiments, the turbulence statistics were measured with LDV using a FiberFlow system (Dantec dynamics, Skovlunde, Denmark) with a 2-velocity component configuration in coincidence mode and 250 mm optics. The LDV measurement volume was traversed in the symmetry plane through the same channel section, which was measured with MRV. The coordinates of the traversing system were calibrated by a simple target placed in the symmetry plane of the channel at a known position. The step size of the traversing system is 0.5 mm in y-direction and six lines in x-direction were measured. Due to the limited accessibility of the laser beams into the flow volume, it was only possible to measure positions that were more than 5 mm away from the walls. At each measurement position, LDV data was sampled over 100 s resulting in 15,000 to 85,000 valid coincidence samples.
Since the LDV measurement volume is substantially smaller in z-direction than the thickness of the MRV measurement slice (6 mm), additional LDV measurements were performed at z = − 3 mm and z = 3 mm to examine the flow variations inside the MRV measurement volume. The variations were not measurable, and therefore, only LDV measurements performed at z = 0, were taken for the comparison with MRV.
The pure acquisition time of all data points, which were used for the comparison with MRV, resulted in more than 15 h. It is interesting to note that this time is more than double the acquisition time of MRV although much less RST data points were acquired (1,400 in LDV versus 90,000 in MRV).

LDV data processing
The LDV processor produced two mean velocity components u x , u y , and three RST components In addition, two higher statistical moments, the skew and kurtosis, were derived from the velocity distribution. The velocity samples were weighted by their transit times to remove the bias between slow-and fast-moving particles in the flow field. The calculation of the velocity moments was done in BSA Flow Software (Dantec dynamics, Skovlunde, Denmark). Additional computations were conducted in Matlab 2018a (The Mathworks, Natick, USA).

LES study
A wall-resolved LES was conducted in the same periodic hill geometry to provide a second reference data set for the validation of the MRV technique. The LES was performed with the commercial CFD analysis code, Star-CCM + 2019R2 (Siemens PLM Software, Plano, USA). The same criteria for generating a wall-resolved LES mesh as described in Addad et al. (2008) were followed. Furthermore, the Wall-Adapting Local Eddy-viscosity (WALE) sub-grid scale (SGS) model was used (Addad et al. 2008). The advantage of the WALE SGS model is that it does not require an ad-hoc near wall damping function but directly handles the wall treatment in the LES domain.
Only the span between two hill crests was considered for the computational domain. Streamwise periodic boundary conditions were used to simulate an infinitely long channel. All other boundaries were assumed smooth and the no-slip boundary condition was applied at all walls. The Reynolds number had a fixed value, which required adjustment of the mean pressure gradient over time to keep the mass flux constant. These well-accepted practices are described in more detail in Breuer et al. (2009).
The resulting mesh had 25.6 million cells and 15 nearwall prism layers to ensure a good near wall resolution. The boundary conditions for the CFD analysis and associated flow properties were determined from the experimental test condition and are listed in Table 1. The Kolmogorov timescale of the problem, estimated from a priori RANS simulation, is 5.21 ms. This scale was used as a guidance for the numerical time step, which was set to 1.00 ms. More information about the LES parameters are given in Table 3.
After performing the LES for a physical time of 1 s, the turbulent flow was considered developed and the velocity data was sampled over a period of 3.5 s, which is equivalent to 12.4 convective times. See Sect. 3.1 for the explanation of the convective timescale. It is assumed that with more than 12 convection times, the qualitative distribution of turbulence in the flow has become clear enough. However, a much longer sampling would be required to achieve fully statistically converged results. Due to limited computational resources, the LES results were taken as they are and used for a qualitative full-field comparison to the MRV data. Note that the sampling time is orders of magnitude smaller than the minimum sampling time of the MRV data (2230 s) and LDV data (100 s).
Both the simulation and the experiment were compared in an assumed fully periodic state. Small differences between the two data sets may occur due to the surface roughness in the experiment and other manufacturing inaccuracies. It is assumed that all errors in the boundary conditions are small and do not affect the qualitative comparison between MRV and LES.

LES data processing
Since the MRV data produced only results in the symmetry plane, the same plane was used for the post-processing of the LES data. The moments of the velocity field, i.e., mean velocity vector and RST, were computed directly from velocity samples that were obtained from all grid points in the symmetry plane for all time steps. It is worth mentioning that the amount of TKE modelled in the SGS of the LES is less than 1% of the calculated TKE, which means that most of the turbulent eddies are resolved in terms of temporal and spatial scales. All processing steps were conducted in ParaView 5.0 (Kitware Inc., Clifton Park, NY, USA).

Results
The following paragraphs present the 2D6C RST data in the periodic hill channel which was reconstructed from the MRV magnitude images. Two MRV reconstruction methods are compared (Sect. 4.1). The final MRV data is then compared to the LES (Sect. 4.2) and LDV data (Sect. 4.3). Finally, the errors in the MRV data are analyzed (Sect. 4.4). Figure 2 illustrates the reconstruction procedure applied in this study. Examples of the raw magnitude images are shown in Fig. 2 A. The depicted magnitude data corresponds to all 12 m enc 1 values for the third row of the encoding matrix in Eq. (2). The signal attenuation strongly correlates with m enc 1 . Also, it can be seen that the signal attenuation is highly non-uniform: The magnitude values in the shear layer downstream of the hill start attenuating with low m enc 1 and fall quickly below the noise floor for larger encoding values. In comparison, signal attenuation in the free stream region is only visible for the highest m enc 1 values. For the demonstrated case, it is difficult to select a single m enc 1 value that leads to sufficient signal attenuation in all flow regions, especially if the turbulence range is not known a priori.

Comparison of MRV reconstruction methods
This finding is supported quantitatively by Fig. 2b. The graph shows the magnitude values for the two marked positions in Fig. 2a, one in the shear layer (red marker) and one in the freestream region (blue marker). Each point corresponds to the mean of the magnitude from all 256 acquisitions. The depicted confidence intervals equal the standard deviation of all acquisitions. Overall, the standard deviation is lower than 0.03 and the Gaussian behavior is visible in the magnitude distribution.
Furthermore, Fig. 2b shows the three-dimensional Gaussian fit over all six encoding directions as described in Sect. 3.1. Note, that here, only the one-dimensional projection on the second encoding direction is shown. It can be seen that the Gaussian fit matches well with all data points. As described in Sect. 3.2.1, all magnitude data points that are below 10% of the normalized magnitude were excluded from the analysis to avoid data points that contain a magnitude bias due to rectified noise. Figure 2c and d compares the Gaussian fit with the case that the three-dimensional Gaussian would be calculated with a single m enc 1 value. For demonstration, the Gaussian functions are calculated with the data for m enc 1 = 15mT m ⋅ ms 2 and m enc 1 = 45mT m ⋅ ms 2 . For ease of interpretation, the data are shown on a logarithmic scale over the square of m enc 1 , since the logarithm of a Gaussian function is a linear function. It can be seen that the slope of the lines depends strongly on which m enc 1 value is used. In the case of high turbulence (orange marker), the Gaussian function corresponding to m enc 1 = 45mT m ⋅ ms 2 leads to substantial deviations, which results in up to sevenfold underestimation of the RST components compared to the Gaussian fit to all 12 m enc 1 values In summary, the three-dimensional Gaussian fit over all 12 m enc 1 values and all six encoding directions results in more

Qualitative comparison between MRV and LES data
Prior to the comparison between the CFD and MRV results, the periodicity of the flow is verified with the MRV data. Figure 3b shows the profiles of the velocity magnitude and Fig. 3c shows the turbulent 2 at the position of the 12th and 13th hill crest. Taking into account the 1-pixel shift in y-direction between the successive hills, which is caused by a slight misalignment of the image axes, the corresponding Pearson linear correlation coefficient for the velocity magnitude and turbulent kinetic energy is 0.97 and 0.98. This indicates that the periodic flow is fully developed at the 12th hill and that the periodic boundary conditions of the LES are indeed valid. Figure 4 shows the comparison between the mean velocity results from MRV and LES for the symmetry plane. While the two data sets show an excellent agreement in the u x component, qualitative deviations are visible in the u y and u z data, which appear simply because of the smaller depicted value range. Figure 5 shows the comparison of all six Reynolds stresses. The LES obtained contours appear noisier because the velocity moments have not been fully statistically converged. The turbulent structures are still visible in the averaged data. The MRV data is remarkably smooth, which underlines the high stability and dynamic range of the three-dimensional Gaussian fit over multiple m enc 1 values, as explained in Sect. 3.2.1. Overall, the two data sets agree well and there are no major qualitative differences visible. The same turbulent features appear in both data sets, for example, the turbulent region downstream of the hill and the strongly anisotropic turbulence at the slope upstream of the hill's crest (visible in u ′ z u ′ z ). As expected, the components u ′ x u ′ z , u ′ y u ′ z as well as the mean velocity u z are close to zero everywhere in the flow because of symmetry.
Small qualitative differences between the two data sets can be observed in the shear layer close to the hill. The MRV data in Fig. 5 D shows a small region with positive values close to the hill crest, which does not appear in the LES data. In a qualitative comparison to similar experiments (Rapp and Manhart 2011), it can be concluded that this is a measurement error in the MRV data. The LES results show the correct behavior.

Quantitative comparison between MRV
and LDV data Figure 6 shows the quantitative comparisons between the available LDV data and the corresponding MRV data extracted at the same position in the flow field. The position of the extracted lines is shown at the top of Fig. 6. As described in Sect. 3.2, MRV captures velocity data directly near the wall, whereas LDV requires a certain minimum distance due to the restricted optical access. In the regions where both data are available, the mean velocity results are in very good agreement with the exception of a few regions near the crest of the hill. However, these deviations should not be given much importance as the maximum deviation (0.04 m/s) is only 3% of the axial bulk velocity at the crest of the hills. Larger deviations are visible in the Reynolds Stress data. As the most striking feature, it can be seen that the u ′ x u ′ x component shows relatively large quantitative deviations in the shear layer downstream of the hill: the peak stresses in the MRV data are up to 0.06 m 2 /s 2 higher than in the LDV data. All positions outside this region show a much higher agreement. A similar trend is observed in the Reynolds shear stress component u ′ x u ′ y . In comparison, the u ′ y u ′ y component shows a substantially higher quantitative agreement in all flow regions. The values of the peak stresses coincide within 0.014 m 2 /s 2 . The reasons for the seemingly directiondependent behavior are analyzed in the next section.

Error contributions in the MRV data
The previous sections showed a good qualitative agreement between the MRV and the reference data sets. However, some quantitative deviations appeared locally, which are here analyzed in detail.

Random errors
In this study, the MRV data was heavily averaged such that it can be assumed that the remaining random errors are negligible. The clean contours in Figs. 4 and 5 and the smooth profiles in Fig. 6 suggest this assumption. All deviations observed in this study between MRV and LDV results likely come from systematic error sources. Note that quantifying the random error, hence the measurement uncertainty, is a complicated task in this work because of the non-linear Gaussian fit of the image magnitude data.

Temporal filtering caused by the finite encoding process
As pointed out by Gao and Gore (1991), Gatenby and Gore (1994) and more recently by Elkins et al. (2009), there is a lower limit of turbulent timescales that can be measured bias-free using the reconstruction technique outlined in Sect. 2. Strictly speaking, the relation in Eq.
(1) is only valid if the Lagrangian integral timescale of the turbulent flow is much larger than the duration of the velocity encoding process, which is the length of the bipolar velocity encoding gradient (see Table 2). An estimate for the Lagrangian integral timescale is here obtained from the LDV velocity statistics. Note that LDV obtains the velocity fluctuation from a fixed position in space, hence the Eulerian frame of reference. Nonetheless, the integral timescales obtained from such measurements can serve as a conservative estimate of the lower limit of the Lagrangian integral timescale (Shlien and Corrsin, 1974).
As a well-known relation, the integral timescale equals the integral of the temporal autocorrelation function divided by the velocity variance. Since the LDV data is measured with random samples, a re-sampling of the data is necessary, which is here conducted with arrival-time quantization (Damaschke et al. 2018). Figure 7 depicts the derived integral timescale for a stream-wise position shortly downstream of the hill, at which relatively short timescales can be expected. The smallest integral timescales occur in the shear layer, where the estimated values are approximately twofold larger than the velocity encoding time of the MRV sequence. Since the integral timescale is derived in a very conservative way, the true ratio is expected to be larger. Such results support the assumptions made in Eq. (1) is valid throughout most of the flow except for regions closer to the hill, where the shear layer is much thinner and the ratio of timescales may approach values closer to 1. For this reason, it seems unlikely that the filtering effect of the finite encoding process is responsible for the large deviations in the u ′ x u ′ x component.

Spatial filtering and bias caused by the finite voxel size
The formula in Eq. (1) is derived assuming homogeneous flow conditions in each voxel. As Kuethe and Gao (1995) pointed out, this assumption can produce significant errors in flow regions with relatively small spatial scales. In particular, the variation of the mean velocity within a voxel leads to additional intra-voxel phase dispersion that results in overestimation of the turbulence quantities. For the presented experiment, the effect of the voxel size is considered small for most of the flow region since the geometry of the flow experiment is large compared to the voxel size in x and y-direction. Also, the LDV results indicate that the variations of the Reynolds Stresses and the mean velocities within the slice along the z-direction are sufficiently small. Errors caused by intra-voxel flow variations are assumed negligible for the streamwise flow positions shown in Fig. 6. More significant errors presumable occur closer to the hill crest where the thickness of the shear layer is smaller than the voxel size. The comparison to LES data in Fig. 5d showed a qualitative deviation in this region, which might be related to inhomogeneous flow conditions within these voxels.

Bias caused by higher orders of motion
The MRV method used in this study is based on a GRE sequence, which can be considered as the standard technique for fluid mechanics experiments with MRV. The weak point of this sequence design is the readout gradient, which is used for the so-called frequency encoding of the space. In comparison to purely phase-encoded sequences known as single point imaging (SPI) methods (Emid and Creyghton, 1985;Mastikhin 1999;Balcom, 1996;Bruschewski et al. 2019), the readout gradient makes the image phase more sensitive to higher orders of motion, which is difficult to control. Particularly noteworthy is the sensitivity to acceleration, which can be related to the second moment of the gradient waveform m 2 . Similar to the velocity dephasing related to m enc 1 , signal attenuation also occurs because of acceleration dephasing and non-zero m 2 . The same effect applies to all other orders of motion. For example, a strong signal attenuation related to jerk dephasing and the third gradient moment m 3 could be observed in John et al. (2020).
In the present sequence, m 2 and higher gradient moments are particularly large for the frequency-encoded x-direction. Note that the calculated m 2 in the x-direction is 1.8-fold larger than in the phase-encoded y-direction and this ratio can be larger for higher moments. Also, the acceleration Fig. 7 Comparison of the velocity encoding time of the MRV measurement (blue) and the integral timescale obtained from LDV (red) at a position inside the shear layer in x-direction is relatively high. These considerations agree with the observed deviation between the MRV and LDV Reynolds Stress data: All components which depend on the fluid motion in the x-direction ( u ′ x u ′ x and u ′ x u ′ y ) show large deviations, while u ′ y u ′ y shows much better agreement. Therefore, it can be assumed that dephasing effects due to higher orders of motions are a major error contribution in the presented MRV data.
While this bias is a significant source of error in the presented work, its influence may well be lower in other implementations, especially if SPI methods are used.

Bias caused by the non-Gaussian turbulence
As mentioned in Sect. 2, the underlying assumption of the MRV reconstruction technique used in this study is that the probability of the random velocity fluctuations must follow a Gaussian distribution. This assumption is often made, but the kind of induced error and its extent is not clear.
Since a Gaussian fit was performed in this work for more than two magnitude data points along each direction, it is possible to estimate the residual error of the Gaussian fit. A qualitative measure is obtained from the residual sum of squares of the three-dimensional fit, which is normalized by the number of fitted m enc 1 encoding points and the maximum value in the flow field. The results are shown in Fig. 8a. It can be seen that the distribution of the residual sum of squares agrees well with the local deviations between MRV and LDV: The highest residual sum of squares occur at positions with the highest deviations in the Reynolds Stress data (see Fig. 6). Figure 8b and c shows how the MRV data would perform if the underlying data was indeed Gaussian distributed. For that reason, the LDV data is processed with the same assumptions as the MRV data. Hence, the LDV velocity samples were first fitted onto a Gaussian distribution instead of calculating the velocity variance directly. As seen in Fig. 8c, the deviation between Gaussian-fitted LDV and MRV data is almost non-existent for u ′ y u ′ y , which is a remarkable result and a strong indication that the Gaussian fit is responsible for the errors in the MRV data. In the case of u ′ x u ′ x , the Gaussianfitted LDV data is only slightly closer to the MRV data than the true LDV data, see Fig. 8b. In this case, other sources of errors dominate, and the most dominant one is presumably related to higher orders of motion in combination with the applied sequence as described in Sect. 4.4.4

Discussion
In theory, the most accurate MRV method for turbulence quantification is Fourier velocity encoding, which provides the velocity spectrum within a voxel without placing assumptions on the shape of the velocity distribution (Moran 1982;Callaghan et al. 1988;Li et al. 1994;Newling et al. 2004). However, this method typically results in long measurement times when a large number of m enc 1 encoding points is used. A less accurate but much faster method calculates the Reynolds stresses from a single m enc 1 encoding value. The approach presented in this study resembles an intermediate step between single m enc 1 encoding and Fourier velocity encoding, hence, a trade-off between measurement accuracy and measurement time.
The presented method and the single m enc 1 encoding method share the same principles. Both techniques are based on the relation in Eq. (1), which describes the velocity spectrum as a Gaussian distribution. With this model, Fig. 8 Influence of the Gaussian assumption on the MRV data reconstruction. A: residual sum of squares of the Gaussian fit in the MRV data, normalized by the number of fitted m enc 1 encoding points and the maximum value. b, c Visualization of the error caused by the Gaussian assumption by comparing the MRV data (blue), the LDV data (red), and the LDV results obtained from Gaussian-fitted raw data (green) the Reynolds stress data can be obtained from the attenuated signal in each voxel, which reduces the measurement time substantially compared to Fourier velocity encoding. However, the two reconstruction steps, i.e., Gaussian assumption and correlation with signal attenuation, are responsible for various types of measurement errors. Other errors arise from the parameters and design of the pulse sequence. In principle, all the identified error contributions can be resolved, or at least substantially reduced, by Fourier velocity imaging in combination with a purelyphase encoded pulse sequence (Callaghan et al. 1988;Newling et al. 2004). However, such a technique is often not applicable due to the very long measurement time.

Achieved improvements using the presented method
The presented method improves the dynamic range and precision of the single m enc 1 method. Based on the ICOSA6 encoding scheme used in Haraldsson et al. (2018), the encoding process applied here is comprised of six encoding directions and multiple encoding values. For each voxel, all six RST components are derived from a single three-dimensional Gaussian fit. Such measurements take longer than Reynolds stress measurements with a single m enc 1 value, but they have several advantages: • The dynamic range is increased in comparison to the RST reconstruction from a single m enc 1 value (see Fig. 2 B + C). • The quality of the Gaussian data fit may provide information on non-Gaussian turbulence (see Fig. 8a), which could not be determined from a single m enc 1 value. • The measurement efficiency (data and precision per time) is higher than the RST reconstruction from individual encoding directions because information is shared between encoding directions in the ICOS6 method (Zwart and Pipe 2013).
With the current method, Reynolds stress data can be obtained in highly complex flow geometries with low uncertainty in an acceptable time. It should be noted that in this study the Reynolds stresses in a strongly non-uniform flow were reconstructed from only 12 m enc 1 values. A Fourier velocity encoded measurement would require a much higher number of encoding points to achieve the same velocity resolution.

Recommendations for further developments
The identification of the two major error contributions, higher orders of motion dephasing and non-Gaussian turbulence, opens up new research possibilities. In addition, several other minor issues and potentials were identified in this study. Future research may focus on these points: • Improved statistical model: The effect of a non-Gaussian distribution is a known problem in Diffusion Tensor imaging. A common method is to perform a onestep Taylor series expansion of the exponent in Eq.
(1), which adds another parameter to the curve fit and accounts for a non-zero excess kurtosis (Jensen et al. 2005). This method could be applied here to reduce the bias by non-Gaussian turbulence. • Noise compensation techniques. The rectified noise in the voxels leads to an overestimation of the signal intensity, which becomes significant if the signal is close or smaller than the noise floor. For this reason, these data points were excluded from the data fit in this study (see Fig. 2 B + D). Using these points would increase the dynamic range of the measurement. • Quantifying the statistical error in the final data. The uncertainty estimation must consider the image noise and the uncertainty of the data fit.
Furthermore, the findings of this work are linked to the highly anisotropic and inhomogeneous turbulence in the investigated setup. The error contributions may differ in other setups and therefore further studies are needed to investigate in how far these findings can be generalized.

Conclusion
This study presented MRV Reynolds Stress measurements in a periodic hill channel with a hill Reynolds number of Re = 29,500. A novel measurement routine was applied, in which MRV data sets are measured with multiple encoding points and encoding directions. A three-dimensional Gaussian fit over all data sets yields all six Reynolds Stresses. The MRV results were compared to a wall-resolved LES. Furthermore, LDV turbulence measurements conducted in the same channel served as the ground truth.
It was shown that the MRV Reynolds Stress data have excellent precision and agree qualitatively with the reference data sets. However, there are apparent systematic deviations. One of the most prominent error contributions is the signal attenuation caused by higher orders of motion, which depends on the setup and on the type and parameters of the MRV method. A fundamental error was also identified in the common assumption that the turbulence is Gaussian distributed. With the presented reconstruction technique, the MRV data are fitted to a statistical model, and depending on the examined flow setup, the Gaussian model can lead to considerable errors. Possible ways how to reduce these errors were presented.
In summary, the presented MRV method enables rapid turbulence measurements in complex internal flows. Compared to other experimental techniques such as HWA, LDV, PIV, and PTV, measurements can be carried out much more easily as no optical or physical access is required. However, there are several sources of error that are specific to MRV and further development work is required.