MRV challenge 2: phase locked turbulent measurements in a roughness array

Magnetic resonance velocimetry (MRV) provides capabilities to measure three velocity components across three-dimensional fields of view without the need for optical accessibility. Predominant usage of the technique in engineering applications centers on steady flows and the discussion of a single time-averaged set of data. In this second MRV Challenge comparison, phase-locked MRV measurements of a pulsatile flow are conducted by five teams from across the globe. A geometry is explored which consists of turbulent flow past six identical cubic elements placed in the center of a square water channel, with partial elements along each side of periodically varying heights. A pulsatile injection flow through the floor between the second and third cubic roughness elements creates a three-dimensional jet that interacts with the complex cube wakes. The details of each scan sequence are summarized for the teams. Results are compared for three different time windows—one pre-injection, a second while the injection flow accelerates, and a third during a quasi-steady condition with the injector fully on. Finally, the influence of the temporal resolution selection for the phase locked MRV is discussed. The results are remarkably similar despite the complex flow configuration selection and showcase a relatively underutilized capability to obtain time-dependent, volumetric data for periodic flows through three-dimensional geometries using MRI. Recommendations for best practices are also provided.


Introduction
In the 2019 MRV Challenge ), phase contrast scans representing a single time averaged velocity field for turbulent flow through a U-Bend geometry were conducted by several international teams. Each team used different magnetic resonance imaging (MRI) systems, sequences, and software parameters to measure three-component (3C) velocity fields in a three-dimensional (3D) domain for a single geometry operated at the same flow conditions. Both commercial and custom pulse sequences were used, and k-space data were reconstructed and filtered by groups using a variety of techniques. The methods and hardware each team leveraged were summarized, and the data from each team were compared in several ways that highlighted the flow three-dimensionality with excellent agreement between the different groups in a geometry that would have made other optical and probe-based measurements difficult.
This second MRV Challenge again considers a turbulent flow field and geometry common to the participant groups but adds the requirement for triggering the MR system to 1 3 28 Page 2 of 19 scan at intervals throughout a periodically pulsed flow. Triggered or gated MRV scans produce phase-averaged, 3D, 3C velocity fields for multiple time intervals or phases in the periodic cycle of the flow (Thompson and McVeigh 2004;Elkins and Alley 2007;Töger et al. 2016;Wymer et al. 2020). In these gated scans, a complete set of k-space data is constructed for each phase of the cycle with subsets of the full k-space data acquired over successive cycles. Each phase-average, therefore, represents the velocity field during a specific phase interval of the cycle averaged over 100-1000 s of cycles. Gated phase contrast scans are commonly used with electro-cardiogram (ECG) hardware to measure cardiovascular flows in biologic systems (Pelc et al. 1994; Thompson and McVeigh 2004;Wymer et al. 2020). In this study, the gating signal was supplied as an analog signal from a data acquisition system (DAQ). The number of phases acquired depended on the MR sequence used and the compromise between temporal resolution, spatial resolution, total scan time, and signal-to-noise ratio. Because the flow is turbulent, complete scans were repeated several times and averaged to reduce measurement uncertainty on the phase-averaged velocity. The final data provide temporal information about the periodic flow from a pulsed injection and, as such, highlight the use of gated MRV for flows with a strong periodic component. The objective of the challenge is to benchmark consistency across independent research groups, to assess measurement uncertainty, and to suggest a set of best practices by identifying commonalities and differences in scan protocols in relation to the results.
There are plenty of applications where time-dependent flow phenomena are of great interest. Within the context of applications of MRI, these include cardiovascular and respiratory flows (Bammer et al. 2007;van Ooij et al. 2016;Banko et al. 2016;Schiavone et al. 2021;Demir et al. 2022;Li et al. 2022), microfluid oscillators (Badilita et al. 2010;Wassermann et al. 2013), gas turbine cooling (Holloway et al. 2002;Joo and Durbin 2009;Talnikar et al. 2017;Fan et al. 2018), and transient pollutant dispersion within urban canopies (Pullen et al. 2005;Hertwig et al. 2019;Homan et al. 2021). The selection of an appropriate geometry for use in a global challenge included consideration of experimental apparatus size so that logistical challenges of shipping and setup could be minimized, engineering relevance to include turbulent flow applications that can leverage the advantages of MRI techniques, and a desire to include two interacting flows-a steady main flow coupled with a pulsing injection flow that operated on a consistent duty cycle. The flow geometry chosen had to warrant an appropriately turbulent flow with adequate levels of complexity that also met the desired objectives of the challenge: the development of both a series of best practices, as well as quantifying experimental uncertainty for the temporal investigation. The vast majority of gated phase contrast work has been done on physiological flows. However, such flows often contain smooth wall flow separation or laminar-to-turbulence transition (since the cycles often contain periods of near zero flow rate), which could pose additional unwanted challenges for reproducibility across groups. Therefore, physiological flows were not considered in the design phase of the challenge.
Flow over sections with regularized roughness has been well studied for some time and poses a difficult and relevant engineering application. Contexts range from flow in urban areas (Cheng and Castro 2002;Inagaki and Kanda 2008;Heist et al. 2009;Herpin et al. 2018;Wang and Benson 2020), to fundamental wall roughness studies (Roberts 1994;Coceal et al. 2006Coceal et al. , 2007Flack and Schultz 2014), to assessing the impacts of additive manufacturing on device performance (Snyder and Thole 2020;Obilanade et al. 2021). In all cases, experimental measurements can be hindered by the complexity of the roughness elements, because these elements may limit the domain where measurements can be made, even though understanding the flow within the roughness layer can be important. MR-based measurement diagnostics are therefore well suited to such flows since they have the added advantage of measuring the flow between and over appropriately sized roughness elements and provides simultaneous information on their location and size.
In this second iteration of the MRV Challenge, an investigation into the turbulent flow through a regularized cube array is conducted that includes a pulsatile surface injection in the wake of a central roughness element. Five teams from around the world conducted gated phase contrast scans using a variety of MRI systems and post-processing techniques, while sharing the flow apparatus and pulsatile pump. Each team produced 3D-3C, phase-averaged data sets which consist of 3D-3C flow fields at regularly spaced phase intervals throughout the pulsing cycle, where each phase contains several million measurement voxels. The flows are compared for several phase intervals including injector off, accelerating injection, and injector on (quasi steady). The purpose of the challenge includes showcasing the sensitivities associated with 3D experimental comparisons in a time-dependent flow as well as beginning to establish some standard techniques for conducting these investigations. The work may encourage increased adoption of gated phase contrast scans using MR systems in engineering application-oriented flow fields. Establishing the method efficacy includes discussion on uncertainty for both spatial and temporal components as well as traditional error sources associated with MR techniques. The intended audiences for this paper are scientists and engineers with interest in fluid mechanics and who can consult with MRI specialists to assess or conduct measurements of turbulent periodic or steady flow fields within the limitations of the technique.

Experimental setup and conditions
The geometry consisted of an upstream flow conditioning section, a test section, and an injector as shown in Fig. 1a. This apparatus was shared by each research group. The flow conditioning section is the same as was used in the first MRV challenge (see Benson et al. 2020, for details). A series of honeycomb elements, grids, and a 4:1 area ratio contraction eliminates secondary flows and creates a nearly uniform inlet flow to the test section. Extensive flow conditioning was important since each team separately determined the flow path to and from the main apparatus based on requirements of their respective MRI facilities. Each component was assembled precisely using alignment pins to ensure repeatability between the teams.
The test section is channel with a square cross section of height H = 25 mm. It begins with a 1 mm square boundary layer trip along all four channel walls and positioned 37.5 mm downstream of the end of the contraction. Six cubic roughness elements are placed along the bottom wall and on the channel center plane, with partial roughness elements on each sidewall as depicted in Fig. 1b, c. The roughness elements begin 41.5 mm downstream of the boundary layer trip. Each cubic element has a height of 7.5 mm and is spaced 15 mm apart. The sidewall elements have the same streamwise spacing, extend 1.75 mm into the flow, and alternate between 7.5 and 12.5 mm tall. Between the 2nd and 3rd element, a hole in the channel floor allows for the injection of a secondary stream. The hole has a 5 mm wide square cross section and is centered between the cubes. The injector is affixed to the channel floor from underneath. It has a connection fitting that runs parallel with the test section, contracts and turns 90 degrees into a vertical orientation as depicted in Fig. 1c and runs straight for 9 injector tube widths before it intersects the test section floor. The coordinate system throughout the results is located on the spanwise center of the downstream side of the boundary layer trip on the test section floor (visible in Fig. 1b). The streamwise coordinate is x, wall normal is y, and spanwise is z. The distance from the origin to the upstream side of the first building is 41.5 mm (1.66 H).
Flow to the main and injection flows were supplied by separate pumps and recirculated from a main reservoir. Each team independently supplied a main flow pump, the reservoir, main flow inlet and outlet tubing, and additional metering components such as flow meters and valves. The set-ups particular to each team are described in the following sections. The pulsatile injection pump along with the tubing from the pump to the injector was shipped with the main apparatus to improve consistency in the injection flow supplied by each group. The pump was a Cole-Parmer ISM446B gear driven motor with Micropump 81,281 pump head. Fig. 1 Overview of the flow conditioning and test section geometry. a Flow conditioning and test section isometric cut view highlighting key components, the main flow, and the injection flow. b Isometric cut view of the test section indicating the boundary layer (B.L.) trip, the sidewall and central cubic roughness elements, the location of the injection tube, and the coordinate system. c Dimensions in millimeters of the test section geometry in section view on the channel center plane The working fluid for all teams was water with a dilute solution of copper sulfate. The main flow rate was 22.2 L per minute, corresponding to a bulk velocity of 0.6 m/s. This produces a main flow Reynolds number of Re = 15,000 based on the channel height, H, the bulk velocity, and the kinematic viscosity of water at 21 degrees Celsius. This flow rate also ensures that the boundary layers produced by the trip are fully turbulent, as described in Benson et al. (2020). The injection flow rate was pulsatile, ranging from 0 L per minute when fully off to 1 L per minute when fully on. In the fully on condition, the bulk velocity through the injector is 0.67 m/s. This gives a blowing ratio of about 1.12, defined as the ratio between the bulk velocities in the injection and main flows. The injector Reynolds number based on its bulk velocity and width is 3400.
The injector flow rate changes from zero to 1 L/min. Since every team used the same apparatus, the flow should be the same aside from the slight differences in the beginning and end of the waveforms due to differences in the impedance of the flow setup. These differences likely overwhelm any small variations that occur because of slightly different accelerations through the transitional regime of the flow.
The injection waveform was generated by filtering a square wave of 1.5 s period and 45% duty cycle with a triangular filter kernel of 0.45 s width. This resulted in a nominal injection time of 1.125 s, with the injector fully off during the remaining portion of the 1.5 s period. The time between the mid-points of the rising and falling sides still represents 45% of the period. Filtering the rise and fall of the underlying square waveform was important for reducing overshoot and undershoot in the actual injection flow due to fluid inertia in the injection line. This waveform was sent to the pump using an analog voltage signal. Importantly, while each team used their own pump and tubing for the mainstream flow components, the injector pump, excitation waveform, and tubing from the pump to the injector was shared between groups to remove as many differences as possible.
The gating as described earlier uses the MRI system's ECG trigger and can be achieved using either prospective or retrospective techniques. Each team used a DAQ system to send a rapid TTL pulse to the ECG trigger that was synchronized with the analog voltage waveform used to control the pump. The alignment of the trigger signal within the injection cycle was determined separately by each team to allow for timing considerations with the different MRI pulse sequences.

MRV details by research team
A complete review of MRV fundamentals is outside the scope of this paper, and the reader is referred to Elkins and Alley (2007) and Wymer et al. (2020) for an introduction to the technique. Some familiarity with the basics of MRI and pulse sequences is assumed for the remainder of this paper. All five teams used gated pulse sequences based on a low flip angle short repetition gradient echo sequence. Rostock and the University of Illinois used FLASH on their Siemens systems. Stanford used the SPGR sequence on a GE scanner, and Hanyang used the T1-FFE sequence on a Philips system. All teams defined the frequency direction along the SI direction in the magnet or streamwise (x-axis) in the flow, and the phase encoding directions were in the vertical, y, and spanwise, z, directions. This section provides additional scan details for each of the five research teams who participated in the challenge. Table 1 summarizes the key equipment and parameters for the experiments. Table 2 describes the velocity encoding strategies of each team.

Hanyang University (Hanyang) scan details
The dataset of the Hanyang University group was obtained using an MRI system located at Ochang campus of the Korea Basic Science Institute in Chungju, Korea. Commercial three-dimensional, three component (3D-3C) flow MRI sequences embedded in the scanner were utilized to obtain the pulsatile velocity data. The sequence used in the experiments consists of velocity encoding using a 6-point phase contrast method and 3D Cartesian orthogonal space encoding. Bipolar gradients were used for the velocity encoding, and the encoding matrix is listed in Table 2a. As illustrated in Fig. 2a, the data for the x component of velocity was sampled first using encoding points 1 and 2 in Table 2a. Point 1 is the reference and point 2 contains a bipolar for the x direction. Next, the data for the y velocity component is acquired, and, finally, the z component is sampled. Flow compensation was applied to gradients in the slice select and frequency directions as well as the gradients in the reference scans that did not have bipolar gradients.
The 5 to 20 mV TTL gating signal generated by a DAQ device (USB-6361 multifunctional I/O device, National Instruments) was transmitted to the MRI scanner through the wireless ECG device (MR400 ECG/Spo2 Transmitters, Invivo Expression) to control the scan in synchronization with the flow cycle. The 1.5 s period was divided into 20 temporal phases to obtain velocity data with a time resolution of 75 ms, and no interpolation between temporal phases was applied.
The main flow was driven using a multi-stage centrifugal pump (PP2-20Y, Hwarang System) and was monitored by a flowmeter (KTM-800, Kometer) located between the test section and the pump to maintain a constant flowrate. The injection flowrate was monitored by a differential pressure flowmeter (L-2LPM-D5V, Alicat) located between the test section and the pump. The voltage waveform was set to have maximum and minimum values of 1.64 V and 0.34 V, and the entire pump system was controlled using the same DAQ device and software with the TTL gating signal for the synchronization of flow and MR scans. The full scans were taken in order of on-off-on. To reduce the background phase error of the spatially corrected data, the process of averaging two flow-on data and subtracting the flow-off data was performed for each voxel. The MRI raw data was reconstructed using the software included with the Philips MRI console. Although parallel imaging methods were not applied to reduce scan time, the zero-interpolation-filling (ZIP) was applied by a factor of 2 in the phase and slice direction. ZIP improves spatial resolution while preserving scan times by zero-filling unmeasured high frequency k-space data, which effectively interpolates in physical space (Zhu et al. 2013). In addition, a house code, based on the algorithm of barrel distortion correction, was applied to correct the spatial distortion caused by the relatively large FOV setting in the reconstructed image. For an accurate velocity mask, the optimal threshold magnitude was determined per slice based on the given cross-sectional area. Finally, the wall-bounded divergence-free smoothing (DFS) filter was applied to reduce the noise in the velocity fields. This filter was modified by the Hanyang University group by adding the no-slip wall boundary condition to the DFS filter developed by Wang et al. (2016).
The uncertainty of the velocity measurement was calculated for each temporal phase and velocity component using the uncertainty equation (Eq. 1) proposed by Bruschewski et al. (2016) where NSA is the number of signal averages, the factor 2 is due to the subtraction of two data sets, and u represents the same velocity component in each data set. The spatial variance is calculated from the difference between two statistically independent velocity data sets comprising all voxels within the region of interest (ROI) of only fluid flow (no walls or surfaces). For the entire temporal phase, the uncertainties did not change significantly and remained within 1.9% of VENC. However, they are slightly higher at the higher injection flow rates (9th to 14th phase) than the other. This phenomenon is probably caused by strong local turbulence and the results are indicative for each scan team.

Stanford University/United States Military Academy (Stanford) scan details
Data were acquired using a system located at the Richard M. Lucas Center for MRI at Stanford University. The triggering system used a DAQ controlled by Labview (National Instruments, Austin, Texas, United States), and the mainstream flow rate was measured by a Transonic Systems TS410 flow module with an ME20PXL probe. The channel was scanned with a 3D-3C gradient echo phase-contrast sequence triggered using the DAQ generated trigger connected to the MRI ECG system. The frequency direction was the primary streamwise direction, x, with a 25.6 cm field of view (FOV), and phase directions were in the coronal direction, corresponding to y, and the sagittal direction, corresponding to z. Bipolar gradients were used for the velocity encoding. As illustrated in Fig. 2b Phase-averaged data were acquired with a temporal resolution of 150 ms using a sequence resulting in 10 phases and a scan time of 20 min and 20 s. These 10 phases were reconstructed in k-space using linear interpolation to give 20 phases for the injection cycle. The ungated flow-off scans were book-ended before and after a flow-on scan, and the averaged flow-off velocity field was subtracted from each phase of the flow-on scan to remove effects from system drift and eddy currents and produce a single phase-averaged set of velocity fields. The 9 sets of these fields were then averaged to produce the final data. A DFS filtering technique was leveraged in post-processing for the final data set (Schiavazzi et al. 2014), which preserves local features while reducing random noise for smoother velocity gradients, helpful in calculating flow quantities where derivatives are used.
The uncertainty was estimated using the technique proposed by Bruschewski et al. (2016) and varies based on temporal phase as a result. The highest uncertainty was estimated to be less than 2% of the VENC and was nearly the same for the vertical and streamwise directions in the region around the jet for the maximum injection phases (9-12).

University of Rostock (Rostock) scan details
Data were acquired using two 3-channel body matrix coils. The mainstream flow was measured by a calibrated SM6020 flowmeter (IFM, Bochum, Germany). The triggering system consisted of a data acquisition unit controlled by Labview (National Instruments, Austin, Texas, United States) using a clock rate of 10 kHz with the digital signal sent every 1500 ms to trigger the scanning. With the same timing, a periodic analog voltage signal was sent to the trigger for the pulsatile pump, which was calibrated in steady-state conditions with a similar flow meter as used for the mainstream.
The minimum and maximum time between two received trigger signals was 1492 ms and 1508 ms. This slight variation was due to the variation in the detection of the trigger by the MRI system and is likely due to the MRI A/D converter and/or the system's gating detection algorithms. The acquisition window size describes the time available for scanning and was set to 1490 ms. This ensured scanning stopped before the subsequent trigger. Each cycle was equally divided into 20 bins with a size of 74.5 ms, which corresponds to the temporal resolution of the phase-locked acquisition. The sequence was designed such that velocity encoding was performed in the inner loop with a loop duration of 4*TR. Bipolar gradients were used for the velocity encoding, and the encoding matrix is shown in Table 2c. As illustrated in Fig. 2b, for each line of k-space, 4 TRs are sampled. Each TR has bipolar gradients on all three axes, and the signs of the bipolars are switched depending on which velocity component is being measured. The 4th TR was the reference, hence, TR 1 and 4 were combined for the x component, 2 and 4 for the y component, etc. Flow compensation was not applied to the gradients in this sequence. The TR was set to 6.2 ms so that the time required for 3 k-space lines fits almost exactly into the bin size (4*3*TR = 74.4 ms). In every 1500 ms cycle, three k-space lines were measured for all four velocity encodes and all 20 injection phases. Each k-space contains 100,080 k-space lines. Therefore, the whole cycle was repeated 3,360 times which resulted in the total acquisition time of 84 min for each flow-on scan.
Relatively low values for the VENCs were used to better resolve the velocities which resulted in some aliasing in the flow-on scans. The aliased velocities were 'unwrapped' using a simple algorithm. Starting from one side of the FOV, the change in velocity from one voxel to the next voxel is compared to the VENC value. The velocity is corrected by adding ± 2*VENC if the change in the velocity value between two voxels is larger than −/ + VENC.
Before and after these scans, a set of three flow-off scans were measured. These scans were not triggered, and the data were averaged and subtracted from each flow-on scan. Due to the lengthy measurement campaign, magnetic field drifts could become significant in the velocity data. It was found that the mean velocities in the flow-off measurements changed by 0.025 m/s between the first and last scans, which were separated by about eight hours. Assuming that this drift is linear, it can be estimated that the velocity error is approximately ± 0.011 m/s over the seven-hour flowon scans. These deviations correspond to 1% of the VENC value and can therefore be considered negligible.
The stochastic measurement uncertainty in the flow-on scans were calculated with Eq.

Seoul National University (SNU) scan details
The data set was acquired at Seoul National University Hospital. The main flow was monitored by a paddlewheel flowmeter. The pulsatile flow pump was operated by a function generator (AFG3152C, Tektronix) using a filtered square waveform with 1.5 s period as an input to the function generator. The flow rate of the pulsatile flow was monitored by a second paddlewheel flowmeter connected to the inlet tube of the pulsatile pump.
A commercial 4-point velocity encoding flow sequence was utilized. A distortion correction filter was applied to prevent gradient warp due to the large FOV. Bipolar gradients were used for the velocity encoding, and the encoding matrix is shown in Table 2d. As illustrated in Fig. 2b, for each line of k-space, 4 TRs are sampled. The first TR is the reference, the second TR contains a bipolar for the x component, and so on. TRs 1 and 2 are combined to construct the x component of velocity, TRs 1 and 3 to construct the y component, and TRs 1 and 4 for the z component. Flow compensation was applied to the gradients in this sequence except the phase encoding gradients. The TTL signal for gating was sent from the function generator to the MRI scanner and synchronized with the waveform signal to the pulsatile pump. The acquisition window size describes the time available for scanning and was set to 1450 ms. Each cycle was equally divided into 22 bins with a size of 65.9 ms, which corresponds to the temporal resolution of the phase-locked acquisition. The sequence was designed such that velocity encoding was performed in the inner loop with a loop duration of 4*TR. The TR was set to 8.21 ms so that the time required for 2 k-space lines fits almost exactly into the bin size (4*2*TR = 65.68 ms). Parallel imaging to reduce measurement time was not applied.
The k-space data were reconstructed using the opensource MATLAB code "mapVBVD" (developed by Philipp Ehses), which reads the Siemens raw data file. Each dataset was then converted to the respective velocity field. The final velocity data were obtained by averaging 2 flow-on scans and subtracting 1 flow-off scan to eliminate background phase error. A volumetric mask was constructed using the signal intensity data to distinguish between fluid and solid regions. The mask was superimposed on the velocity data, forcing spurious data in solid regions to be zero. A DFS filter was applied to the masked velocity data to smooth out noise and reduce error in the velocity field. The DFS filter is based on Wang et al. (2016), and was modified by Hanyang University for wall-bounded flows with no-slip condition at the walls similar to the recommendations in Benson et al. (2020).
Uncertainty in the velocity measurements was evaluated using Eq. (1) by the approach from Bruschewski et al. (2016). The uncertainty, which was estimated from the raw data, is 0.037, 0.032, and 0.016 m/s for the maximum pulsatile flow rate and 0.032, 0.027, and 0.014 m/s for the minimum pulsatile flow rate in the x, y, and z directions, respectively. It is estimated that the uncertainty is higher during jet injection because of the increased turbulence, which is a major source of noise. Nevertheless, the uncertainty values are within 2.5% of the VENC for all directions and phases.

University of Illinois: Urbana Champaign (UIUC) scan details
The MRV experiments were performed at the Beckman Institute at the University of Illinois, Urbana-Champaign. A Blue-White (F1000RB) paddle flow meter calibrated against a Kobold MIM Electromagnetic flow meter was used to control the main flow rate accurately to ± 3.4% of the nominal test flow rate. An Omega FLR1605A differential flow meter upstream of the pulsatile pump, with a 0 to 2 L/min range and a full-scale accuracy of 2%, was used to pre-adjust the pump voltage inputs to ensure a steady flow rate of 0 and 1 L/min with the main flow at nominal conditions. The pulsatile pump was independently controlled by a programmable Keysight 33600A waveform generator that served as a synchronizer providing the damped sinusoidal waveform to the pump and a 5 V TTL trigger signal to the scanner for prospective gating. There was no delay imposed between signals. Two flexible coils, a small 4-channel placed under the test-rig and an 18-channel array wrapped around the top and lateral walls, provided the best signal among the different coil configurations evaluated. The dataset was acquired using a Siemens sequence (fl_pc) with Cartesian sampling, and asymmetric phase encodings with flow compensation for the slice-select and readout gradients. Bipolar gradients were used for the velocity encoding, and the encoding matrix is shown in Table 2d. 4 TRs are sampled. The first TR is the reference, the second TR contains a bipolar for the x component, and so on. TRs 1 and 2 are combined to construct the x component of velocity, TRs 1 and 3 to construct the y component, and TRs 1 and 4 for the z component. Parallel imaging k-space domain based iPAT/GRAPPA with an acceleration factor of 2 in the y-direction allowed for reducing the scan time with negligible impact on the measurements. The acquisition time window was set to 1490 ms, with 5 ms trigger delay at the beginning of each cycle. 31 equally spaced temporal phases were acquired for the same phase encodings. This allowed maximal temporal resolution, limited to the 47.35 ms time (total time for four-point flow encoding and three phase components). No filters or corrections were used during the acquisition or reconstruction.
Non-gated flow-off scans were used to correct for possible eddy currents and system drift over time, none of which were identified over a total testing time of 6 h. Given the absence of drifts, all 8 flow-off scans were averaged for background subtraction, limiting the impact of noise in the flow-off data.
Data postprocessing using in-house MATLAB codes consisting of averaging each temporal phase of the 3 flowon scans following reconstruction, subtracting the averaged flow-off maps, and correcting for biases caused by nonlinearities of the magnetic field gradients used for spatial encoding. The spherical harmonic coefficients for the scanner provided by Siemens were used to model the spatial dependence of the field generated by each gradient coil and used them to correct our data in a two-step process. A 3D image distortion correction was first applied to map the phase-encoded information to the correct spatial location. The correction reproduces the unwrapping method incorporated in the Github Human Connectome Project Pipelines distribution (GitHub). The second step corrects the impact of the gradient nonlinearities on the first moments used to encode the motion, needed for the true magnitude and direction of the velocity data. This phase distortion correction follows the generalized reconstruction technique described by Markl et al. (2003). At last, the DFS filter discussed by Schiavazzi et al. (2014) was applied to each temporal phase to reduce random noise.
The uncertainty in the measurement of each velocity component was estimated from the postprocessed unfiltered data using the approach proposed by Bruschewski et al. (2016), under the assumptions that the noise is constant across the selected region for different realizations, and that its spatial distribution is Gaussian. The voxel-to-voxel difference between two independent velocity maps for the same temporal phase and different scans is used to compute the spatial variance. The temporal phase averaged estimated uncertainties are 2.2%, 2.1% and 2.2% of the VENC values for each direction x, y, z, respectively.

Specific scan timing details
As each team used different values and sequences for reporting echo (TE) and TR, a short description of the general scan sequences is provided to aid investigators interested in duplicating the technique, and to help in overall understanding of the time windows for each since there are fluid dynamics and economic factors associated with overall scan time durations. Figure 2a graphically displays some of the key details associated with the Hanyang University group, who used the 6-point velocity encoding scheme. In their sequence, the velocity is encoded separately for the three directions, by combining a reference excitation with no bipolar gradient applied with a second excitation using a bipolar gradient. This leads to an inner loop time of 2TR. The inner loop time describes the smallest non-divisible time scale in a sequence. Note that all other groups used the 4-point velocity encoding scheme as shown in Fig. 2b, and all velocity information was encoded simultaneously with an inner loop time of 4TR. The inner loop is composed of four alternating bipolar gradients. Combining adjacent bipolar gradient pairs determines the velocity component in each coordinate direction. Therefore, each velocity component is again associated with a 2TR time window. The acquisition of all three velocity components is associated with a 4TR time window. Note that the acquisitions for different velocity components are not collocated in time. Therefore, the temporal resolution of four-point encoding methods is 4TR when considering all velocity components together.
Except for the velocity encoding method, the acquisition process depicted in Fig. 2 is representative of all groups. Specifically, within a single injection period, one or several k-space lines are acquired in each temporal phase. The injection cycle of 1500 ms must be repeated multiple times until all k-space data for all velocity encodings and all injection phases are sampled. Some teams used a relatively small TR so that multiple k-pace lines could be acquired during each encoding segment, which reduces the number of cycles necessary to sample all k-space data. This is the reason for the large deviations in the single scan times in Table 1.
In the case of the Hanyang University study, one k-space line was measured in each segment. To obtain a temporal resolution of 75 ms (20 phases) and given an inner loop time of 2TR (6-point encoding), the TR was adjusted to 37.5 ms. The UIUC team also acquired one line of k-space in each segment but used a slightly different approach. With an inner loop time of 4TR (4-point encoding) and a temporal resolution of 47.35 ms (31 injection phases), the TR was adjusted to 11.8375 ms.
The teams at Stanford University, University of Rostock and Seoul National University leveraged the capability to adjust the number of k-space lines per segment, which means they acquired multiple lines of k-space in each encoding segment. The team at Stanford measured 10 injection phases (reconstructed to 20 injection phases) which resulted in a time resolution of 150 ms. With an inner loop time of 4TR (4-point encoding) and a TR of 5.356 ms, they acquired 7 lines of k-space in each segment. As a result, their scan time was the shortest among all teams. The teams at the University of Rostock and the Seoul National University used a similar approach. They acquired 3 or 2 lines of k-space in each segment, respectively. However, when acquiring multiple k-space lines per segment, the temporal resolution is effectively reduced since the k-space data are spread out over time. For instance, the temporal resolution for one velocity component when acquiring one k-space line per segment is 2TR, while acquiring 3k-space lines with 4-point velocity encoding reduces the resolution to 4TR*3 (or equivalently the length of a segment).
The research grade sequences used by these groups allow for this segmentation, whereas many commercial systems do not have access to this parameter. The cost of multi-segmenting is primarily one of blurring the temporal resolution slightly as the multiple acquisitions do not happen at the precise same time.
Another aspect worth mentioning is the gating technique. With prospective gating, the segments are acquired at exactly the same points in time in the cycle, which reduces blur compared to retrospective gating. However, each trigger signal has a certain amount of uncertainty that must be considered when interpreting the expected gating. Usually, the acquisition window is made slightly shorter than the time between two trigger signals to account for these small deviations. For example, the teams here using prospective gating used acquisition windows between 1450 and 1490 ms (compared to the 1500 ms cycle). As a result, each phase is slightly shorter than the cycle divided by the number of temporal phases and there is a time shift that might accumulate over a cycle. These differences can lead to temporal discrepancies in the data.
With retrospective gating, the k-space lines are acquired with the prescribed timing dependent on the TRs and number of k-space lines in each segment. The gating signal is sampled simultaneously, and the lines of k-space are sorted into their respective phase intervals during reconstruction. If the TRs are not prescribed so that there is an integer number of k-space lines in each phase interval, then the lines start to span across the phase interval boundaries. This is particularly significant for the last phase which depending on the timing may be under sampled. In the case of this experiment, the Stanford TRs were set to match the phase intervals, so this is not a concern. Of course, there is some uncertainty in the timing of the gating signal which may cause some blurring across phases when k-space lines are misplaced during reconstruction. This effect, in general, is negligible relative to the uncertainties described in previous sections.

Results and comparisons
The datasets from each team were linearly interpolated in space onto a common isotropic Cartesian grid of 0.5 mm resolution, which is less than or equal to the resolution used by all teams. Data were also shifted in time to align the start of the pulsatile waveforms, but the native time resolution of each team was preserved. Select times of interest analyzed below are obtained using linear interpolation in time. All other data processing procedures, such as divergence free filtering, were determined and applied by the individual teams as described in the previous sections.

Inter-site variability due to experimental set-up
A challenge for this study was the need to balance reproducibility in the pulsatile flow against the complexity and cost of the set-up in order to transport the experiment between international teams. As discussed above, each team used the same pulsatile pump and inlet tubing for the injector but supplied their own main flow pump and tubing. This resulted in certain flow differences between sites that are ascribed to the selection of these experimental components. In this section, we quantify these flow differences to separate inter-site set-up variability from differences due to MRI procedures. Figure 3 plots the injector flow rate waveforms in liters per minute and normalized by the peak injector flow rate for each team in comparison to the target waveform. Since each team was not able to measure the time-resolved injection flow rate with a flow meter, these results were derived by integrating the MRV data across the injector. This results in uncertainties due to spatial resolution and image segmentation but provides an assessment of relative differences between the teams. MRV data from the Stanford team were compared to time-resolved ultrasonic flow rate measurements and showed good agreement to within about 5% uncertainty, which is easily attributable to MRV resolution. As evident in Fig. 3a, the peak flow rates varied between teams but were typically within 10% of the target waveform during the accelerating and peak injection phases. An exception is the Seoul National University data, which is substantially attenuated. This is likely due to an MRV post-processing step since little data was available in the injector tube for integration, and that waveform should be interpreted qualitatively. We expect that differences in peak flow rate will affect the jet penetration distance and downstream flow evolution, and this should be considered in the interpretation of subsequent results. Significant differences between teams are also observed during the deceleration phases. Specifically, the true waveforms for most teams lag the target waveform. Figure 3b plots the same waveforms normalized by the peak injection flow rate from the MRV data. The normalized plots show good agreement between the acceleration of the jets and the same discrepancy is apparent in the decelerating phases. While the injector pump and tubing were held constant across the teams-the mainstream pump was supplied by the individual teams and ranged in power from 250 to 1100 W between teams. The main flow return tubing length from the outlet of the rig to the reservoirs also varied. Despite efforts to have all teams use similar pumps and main channel tubing and guidelines for tuning the pulsatile pump excitation waveform, the dynamics of the experimental setup differed for all groups due to the main flow component differences. While some groups had less damped waveforms (overshoots in peak and minimum flow rates with quicker transitions) others had damped waveforms (no overshoot but slower deceleration).
In light of the set-up variability, Fig. 3 suggests phases to compare the results from each team. The first depicted phase from each team represents a portion of the cycle where the injector had no flow and would be the baseline comparison condition. This phase is referred to as the Injector Off phase. The rapid change in injector flow preceding the fully on condition presents a condition with local flow acceleration may be of particular interest in the vicinity of the injector and is referred to as the Accelerating phase. This phase is located at 0.4 s into the duty cycle. Finally, the latter phases where the injection is fully on and the 3D phase-locked flow is presumably the steadiest presents a third timeframe to compare. This phase-occurring at 0.82 s into the cycle, will be referred to as the Injector On phase in the subsequent comparisons. Figure 4 depicts the bulk streamwise velocities for each team across the injection cycle. Approximately 20 planes of data from cross sections upstream of the injector and between the first two building elements were used to obtain the bulk velocity from the MRV data at each phase. Specifically, the flow rate was integrated and divided by the known cross-sectional area. The bulk velocity from the MRV data matched the flow meter readings from each team to within about 3% uncertainty, which is the typical uncertainty of the paddle wheel and ultrasonic flow meters used. The Seoul National University data is again an exception, and it is unclear why there is a larger discrepancy between the MRV derived bulk velocity and the independent flow meter reading compared to the other teams. The average velocities in most cases are influenced by the injection pump, with some teams experiencing as much as 5% velocity fluctuations between phases when the injector flow rate is on and off. Differences in the main flow response are due to inter-site variability in the main flow pump power and main flow tubing lengths used, which change the impedance of the system. Because of this behavior, the bulk streamwise velocity for each team obtained from the MRV data between the building pair immediately downstream of the injector will be used to normalize velocity components at the three phases mentioned previously. Specifically, prior to the injection (first phase shown), during the accelerating phase, and with the injector fully on. The latter two timings are depicted with vertical lines in Fig. 4 just as in Fig. 3 previously. This approximately removes some of the differences due to absolute flow rate and provides a fairer comparison between the datasets in order to assess similarities and differences due to the MRI techniques. Figure 5 depicts the normalized velocity magnitude for the first temporal phase from Figs. 3 and 4 for each team at the channel center plane bisecting the cubic elements. Normalization was obtained for this and subsequent velocity magnitude plots using the appropriately timed bulk average velocity shown in Fig. 4 from the MRV data. The geometry is depicted in the gray shaded sections and is imported separately from the data and included to show the edges clearly without the limitations that the experimental resolution would impose. The injector-located centrally past the first solid element in the left portion of the images-has no flow-on and thus the recirculation zones in each element wake are similar. As the boundary layers at the upper wall and above the elements grows downstream-each team presents a roughly similar normalized magnitude increase in the undisturbed mainstream flow as well as the apparent penetration depth into the building wake regions. Similar low velocity regions behind each of the cubic elements are observed in all cases. Since the UIUC team had better temporal resolution than the other teams, the apparent noise in their data is subsequently higher, a feature discussed later.

Accelerating phase comparisons
The selection of an accelerating phase for comparison was conducted by interpolating adjacent phases for each team to a nominal injection time centered at 0.4 s. As is apparent in Fig. 3, this represents a time of rapid injection rate change and presents an important transient interval for flow development especially for the wake regions behind the buildings. Figure 6 depicts this interpolated phase from the MRV data in the same spanwise center plane as in Fig. 5 but includes in-plane velocity vectors which are helpful in visualizing the nature of the flow in the building wakes as well as the behavior of the injected flow. Broadly, the five depictions are generally similar, but there are also some important differences that can be noted. First, the separation at the front edge of the second building varies between teams with little to no separation evident in the SNU data and the largest separation in the Hanyang data. Here-the scan details become important, as the Hanyang data is the best spatially resolved data set of the five teams, and this region is quite small in nature, suggesting that the apparent size may be on the order of the scan resolutions used. The UIUC data has the best temporal resolution as more phases were acquired than the other teams-and may have increased velocity magnitudes on the windward side of the second building as a result. The jet in the UIUC case also appears to be slightly faster than the others-likely due to pressure differences as mentioned earlier. The shear layers between the jet and mainstream flow appear very similar in thickness and velocity vector angles especially in the Hanyang and Rostock images, with subtle differences among the other teams. The jets from each team all bulge in the upstream direction similarly, while the downstream side of the jets have noticeable differences in several cases. The wakes of the 2nd and 3rd building elements in the images show similar magnitudes and recirculation features as well.

Injector on phase comparisons
For teams except SNU and UIUC, the 11th phase from Fig. 3 was used to obtain the quasi-steady injector on condition depicted in Fig. 7. The 12th phase from SNU and the 17th phase from UIUC were used which align nearly exactly with the dashed vertical lines on Figs. 3 and 4. While some of the slight differences in flow features noted during the acceleration phase from Fig. 5 are largely gone; the  injection velocities show variability between all the teams. This is because normalization by the bulk velocity cannot compensate for simultaneous differences in the main flow and injection flow speeds.
These differences are consistent with the injection flow rate differences due to inter-site variability discussed in Sect. 4.1. Despite the jet magnitude differences, there is increased similarity in the size of the separation region above the 2 nd cubic element from Fig. 6. In addition-the region between the injector jet and the building face just downstream also appears much more similar across each team. Similarly, the wake regions after each cubic element are qualitatively in agreement.
To quantify this last comparison, normalized streamwise velocity profiles from each team were extracted where the black vertical line in Fig. 7 indicates the building wake downstream of the injector. Profiles were also extracted from the same location during the injector off phase for comparison. Figure 8 plots the vertical velocity profiles for each team and shows excellent agreement during the injector off phase. The data from most teams lie within the 3% quoted uncertainty band based on the channel bulk velocity, except for the Seoul National University team in the recirculation region near the wall. The injector on phases shows similar agreement despite the significant upward displacement and distortion of the main flow due to the injection flow passing over the cubes. The UIUC profile is slightly more affected by the injection flow, and this is likely due to that data having the largest injection flow rate.
The root-mean-square (RMS) difference across datasets was also computed for a volumetric flow domain. The streamwise extent is the same as shown in Figs. 5, 6 and 7, and it includes the entire channel cross-section. That point-wise variance between teams (i.e., variance at each voxel) was computed for a given phase and averaged over the domain before taking the square root to give a single measure of the RMS difference at that phase. The RMS differences were computed for the injector-off, accelerating, and injector-on phases. The x-velocity (streamwise), y-velocity (vertical), and z-velocity (spanwise) differences were less than 6.5%, 3%, 3% percent of the bulk velocity, respectively, at each phase. Visualization of the RMS difference (not shown) indicated the largest differences at the injection location and near the cuboid surfaces, in agreement with the discussion of Figs. 5, 6 and 7 This supports the consistency of the time-dependent velocity field measured by each team considering the inter-site variability described above.
To better visualize the 3-D nature of the flow and leverage the utility of the three measured velocity components, Fig. 9 depicts the Hanyang University normalized vorticity magnitude of 4.2 for the Injector On phase and includes each of the six centrally located cubic elements within the domain.  Fig. 7 for a the injector off phase and b the injector on phase. The shaded gray region shows the median velocity across the teams with error bands that are 3% of the bulk velocity The vorticity magnitude is calculated from derivatives of each velocity component which are approximated using a second-order central differencing technique built into the software visualization program-TECPLOT™. The normalization was accomplished using the building height and bulk velocity from the same team. The positive z wall has been blanked for clarity while the negative z wall with the partial buildings of alternating heights remains visible. Several key features are evident in Fig. 9. At the extreme lower left of the image-a horseshoe vortex is visible at the base of the leading edge of the first element-common across all groups and the sole element with this readily apparent feature. Vorticity is also evident due to the side-wall boundary layers and along the partial building elements near the visible channel walls. Importantly, the vorticity at the edges of the jet for this phase are clear between the 2nd and 3rd elements which imposes an interesting effect on the downstream building in its wake where the surface is absent the vorticity that occurs at the upstream buildings. This jet effect decays further downstream and the vorticity around the final element closely matches the 2nd building behavior. The ability to measure velocity components at a spatial resolution conducive to computing vorticity directly from the data set is an apparent strength of the measurement technique employed and exists at each measured phase of the flow field.

Temporal resolution considerations
While most teams in the challenge used approximately 20 phases within the injection cycle, one team (UIUC) used 31 phases. This allows for a qualitative discussion on the impact of the number of phases measured for this flow, to include advantages and disadvantages. This section is meant to provide useful interpretation of the tradeoff between statistical convergence of the phase-averaged velocity fields, scan time, and temporal resolution when clinically available sequences have different flexibility in the parameters that can be controlled. For instance, the software may or may not allow for multiple k-space lines to be measured per phase which directly affects a trade-off between scan time and temporal resolution. In order to quantitatively understand the impact of temporal resolution on the flow field reconstruction, future experiments using a single scanner and flow set-up should investigate these differences by systematically varying scan parameters and through comparison to an independent technique, such as phase-locked particle image velocimetry at a region of interest. Figure 10 depicts contours of the vertical velocity component at a cross-plane just downstream of the injector with inplane velocity vectors, using the data from the UIUC team. Figure 10a, c shows contours for the accelerating and injector on (quasi-steady) phases at the full temporal resolution. Figure 10b, d shows data at the same times but temporally filtered by using a moving average of three phases. Therefore, the temporal resolution of the filtered data is approximately three times lower, analogous to measuring 10 phases and reconstructing to 31 using interpolation This is like the procedure used by the Stanford team where multiple k-space lines were acquired per phase, however, temporally filtering data in post-processing is not strictly equivalent to measuring k-space data over a longer time window.
The temporally filtered data during the quasi-steady injector on condition is qualitatively smoother and more symmetric about the center plane than the unfiltered data. In turbulent flows, averaging across multiple phases effectively averages across different turbulent flow structures, which improves statistical convergence of the data. However, averaging may also obscure mean flow transients in time-periodic turbulent flows. For instance, averaging also improves the smoothness and symmetry of the velocity field for the accelerating phase, but differences in the overall velocity magnitude are also evident. The velocities of the filtered data are moderately attenuated relative to the unfiltered data by about 0.1 m/s in the fast, upward moving region of the jet. The temporal resolution of the unfiltered data is 4TR = 47.35 ms since one k-space line was measured per segment by the UIUC team. The filtered data have a temporal resolution of 3*4TR = 142 ms. The acceleration time scale of the jet at this phase is about 225 ms, indicating that the native temporal resolution of the UIUC team well-resolves the acceleration of the jet fluid while the reduced temporal resolution is not much smaller than the acceleration time scale. Note that the temporal resolutions of the other teams range from about 66 to 150 ms (the resolution when acquiring multiple k-space lines per segment is simply the segment length). Therefore, slight differences in velocity fields are expected across teams during the accelerating phases.
The additional tradeoff is that the Stanford scan sequence parameters resulted in a total scan time of 20:20 (min:sec) while the UIUC team had a single scan time of 98:00 (min:sec) which allowed the Stanford team to repeat scans several more times to further improve convergence. Therefore, measuring multiple lines of k-space per phase could be advantageous in slowly varying flows or when scan time constraints are significant. If a sequence only allows for one line of k-space to be measured per phase, then it is optimal to acquire as many phases as possible for a given TR as was the case for the UIUC team.

Conclusions and future work
Five teams across the world conducted MRV measurements in a complex turbulent flow in a square water channel. A freestream Reynolds number of 15,000-based on the channel hydraulic diameter-was similar with the first iteration of this challenge ). The flow passed six identical cubic roughness elements along the center of the channel, with a square injection port located between the second and third elements. Injection happened with a 45% duty cycle to create a pulsatile flow field. Phase averaged three-dimensional flow field measurements were taken by each team for all three velocity components with small temporal windows (phases) ranging from 47 to 150 ms across the cycle. Scan details from each team were assembled for easy reference. The results for the 5 teamswhich included 20-31 data sets (one for each phase) for each team-were compared for three phases of the injection cycle: prior to the injection, the injector turning on, and The left images a, c are the data as measured with the smallest temporal window of the teams; the right have been averaged with adjacent phases to create longer time windows a fully-on injection condition to both understand the flow field and to compare results from each team. The results were largely similar but also showed differences which were assessed to stem from different mainstream pump sizes and tubing as well as the tuning of the injection waveform. These differences as well as differences between the groups' scan techniques (including temporal resolution) impacted the measurements. In order to better compare results from each team, results were normalized using scan specific phase dependent bulk velocities. In addition, the differences related to temporal resolution highlight the tradeoff between temporal resolution and increased scanning time. This second MRV Challenge iteration highlights an aspect of MRV capabilities that has not been as widely leveraged in technological flows-phase locked measurements-which may have broad applicability and extend the use of the technique in flow applications.
Future work in the MRV Challenge will evaluate additional applications of the MR-based techniques beyond velocity measurements, such as the measurement of passive scalar quantities like concentration or temperature. An ability to couple this added capability with increasingly sophisticated velocity measurement techniques will further the applicability of MR-based measurements to complex flows in fluid mechanics and engineering applications. The utility of the measurements in design, research, and assessment of advanced systems is expected to continue to expand as a result.
Authors contributions MB and AB authored the bulk of the text, SS provided Fig. 2, CE interpolated all data onto a common grid for analysis. All authors helped write institutional sections and reviewed the text.

Funding
The overall project was sponsored by the Air Force Research Laboratory (Grant No. A2108-057-021-013385) who supported fabrication costs and some shipping expenses.
Availability of data and materials The data are available upon request from any of the author team provided a host site with sufficient storage is provided with the request.

Conflict of interests
The authors report no competing interests in this work.
Ethical approval N/A.