Jeremías Garay

Hernán Mella

Julio Sotelo

Cristian Cárcamo

Sergio Uribe

Cristóbal Bertoglio
email: c.a.bertoglio@rug.nl

Joaquín Mura

Assessment of 4D Flow MRI's quality by verifying its Navier-Stokes compatibility

4D Flow Magnetic Resonance Imaging (MRI) is the state-of-the-art technique to comprehensively measure the complex spatio-temporal and multidirectional patterns of blood flow. However, it is subject to artifacts such as noise and aliasing, which due to the 3D and dynamic structure is difficult to detect in clinical practice.

In this work, a new mathematical and computational model to determine the quality of 4D Flow MRI is presented. The model is derived by assuming the true velocity satisfies the incompressible Navier-Stokes equations and that can be decomposed by the measurements u meas plus an extra field w. Therefore, a non-linear problem with w as unknown arises, which serves as a measure of data quality.

A stabilized finite element formulation tailored to this problem is proposed and analyzed. Then, extensive numerical examples -using synthetic 4D Flow MRI data * Both C. Bertoglio and J. Mura are joint last authors, listed in alphabetical order.

Introduction

Time-resolved 3D flow magnetic resonance imaging, known as 4D Flow MRI, has shown in the last years potential in assessing cardiovascular diseases since it offers full coverage of the region of interest, therefore allowing for its analysis after the scan [START_REF] Markl | 4d flow mri[END_REF][START_REF] Soulat | 4d flow with mri[END_REF]. Additionally, it allows the computation of several hemodynamic parameters, which can be used as new biomarkers [START_REF] Soulat | 4d flow with mri[END_REF]. However, high-quality 4D Flow in subjects involves long time scans (>20 minutes) even with coarse spatio-temporal resolutions making it challenging for everyday clinical use. In order to accelerate the acquisition time, several strategies have been proposed, such as parallel imaging [START_REF] Pruessmann | Sense: Sensitivity enconding for fast mri[END_REF][START_REF] Griswold | Generalized autocalibrating partially parallel acquisitions (grappa)[END_REF] which accelerates the acquisition by exploiting the sensitivity of multiple receivers, and k-space undersampling [START_REF] Kozerke | Accelerating cardiac cine 3d imaging using k-t blast[END_REF][START_REF] Baltes | Accelerating cine phase-contrast flow measurements using k-t blast and k-t sense[END_REF][START_REF] Gamper | Compressed sensing in dynamic mri[END_REF][START_REF] Ma | Aortic 4d flow mri in 2 minutes using compressed sensing, respiratory controlled adaptive k-space reordering, and inline reconstruction[END_REF] which exploits data redundancies in frequency and time. This scan time reduction comes at the price of reducing the signalto-noise ratio (SNR). SNR also decreases when reducing the image's voxel size. Moreover, the velocity field can only be obtained under a certain predefined range which depends on the magnetic gradient setup, therefore being potentially subject of velocity aliasing. Other artifacts may also appear due to subject's respiration and motion during the scan.

To the best of the author's knowledge, quality control of 4D Flow in clinics is based on calculation of peak/mean flows, mean velocities, flow patterns, and stroke volumes [START_REF] Sabri Koltukluoglu | A robust comparison approach of velocity data between mri and cfd based on divergence-free space projection[END_REF][START_REF] Thekla H Oechtering | Time-resolved 3-dimensional magnetic resonance phase contrast imaging (4d flow mri) reveals altered blood flow patterns in the ascending aorta of patients with valve-sparing aortic root replacement[END_REF][START_REF] Soulat | 4d flow with mri[END_REF]. A more systematic approach is to compute the divergence field of the data: assuming the blood flow is incompressible, jumps in the divergence field may indicate the presence of artifacts. Indeed, the incompressibility assumption has been used for denoising [START_REF] Mura | Enhancing the velocity data from 4d flow mr images by reducing its divergence[END_REF][START_REF] Busch | Reconstruction of divergence-free velocity fields from cine 3d phase-contrast flow measurements[END_REF][START_REF] Sabri Koltukluoglu | A robust comparison approach of velocity data between mri and cfd based on divergence-free space projection[END_REF][START_REF] Ong | Robust 4d flow denoising using divergence-free wavelet transform[END_REF] and as a regularization term during the reconstruction process [START_REF] Santelli | Accelerating 4d flow mri by exploiting vector field divergence regularization[END_REF][START_REF] Bostan | Spatio-temporal regularization of flow-fields[END_REF]. However, an important limitation of the divergence is that any measured velocity will have "infinitely" large divergence compared to its reference value which would be zero. Therefore, to the best of our understanding, only the spatial distribution of the divergence may be used as an error indicator but not in absolute terms for image quality check. Moreover, it will be shown later in the article, measurement artifacts may lead to no important changes in the divergence field.

Therefore, this work introduces an alternative quantitative approach for assessing 4D Flow quality by verifying the compatibility with the linear momentum conservation part of the Navier-Stokes, which when being written appropriately, includes also angular momentum and mass conservation.

The rest of this article is structured as follows. In Section 2, the mathematical model will be introduced, and a numerical method will be developed and analyzed. In Section 3, a set of relevant examples, using synthetic data, will be detailed. Numerical computations of the new model for several types of artifacts and its comparison against the divergence of the data are also shown. Results with phantom and subject's 4D Flow MRI are shown in Section 4. Finally, in Section 5, we discuss potential applications of this metric in the context of reconstruction and processing 4D Flow MRI.

The mathematical model 2.1 The continuous problem

We assume a physical velocity field u, which satisfies the conservation of linear momentum of the incompressible Navier-Stokes equations in the vessel lumen Ω:

ρ ∂u ∂t + ρ u • ∇ u -µ∆u + ∇p = 0 in Ω, ( 1
)
where p is the fluid pressure field, and ρ and µ are the density and dynamic viscosity of the fluid, respectively. We recall that both the convective representation of the advection term and Laplacian representation of the viscous term are written in a simplified form using that ∇ • u = 0. Moreover, the constitutive model (incompressible Newtonian model) represented by the stress terms (viscous plus pressure) is derived by enforcing conservation of angular momentum.

Let us denote u meas the 4D Flow measurement field. We assume that there exist a field w, that satisfies:

u = u meas + w in Ω (2) ∇ • w = 0 in Ω (3) w = 0 on ∂Ω, (4)
with ∂Ω being the whole boundary of Ω. By writing (1) in weak form, and using relations (2)-( 4), we can formulate the following weak problem: Find (w(t),

p(t)) ∈ H 1 0 (Ω) × L 2 0 (Ω) such that Ω ρ ∂w ∂t • v + ρ (u meas + w) • ∇ w • v + ρ w • ∇ u meas • v + µ∇w : ∇v -p∇ • v + q∇ • w = - Ω ρ ∂u meas ∂t • v + ρ u meas • ∇ u meas • v + µ∇u meas : ∇v (5) for all (v, q) ∈ H 1 0 (Ω) × L 2 0 (Ω).
The following remarks are in order.

Remark 1

The left-hand-side of Problem (5) resembles the incompressible Navier-Stokes equation, up to two additional terms. Unfortunately, none of these terms are (semi-)positivedefinite. Therefore, the analysis of the well posedness of this continuous problem (existence, uniqueness, time-stability) becomes challenging. However, those properties can be ensured at the discrete level by including adequate stabilization terms and constraints on the physical constants.

Remark 2 Equation (1) uses the so-called convective form of the advective term. While alternative forms of Equation (5) could be derived starting from other forms for the advection (e.g., conservative), the resulting discrete problem will need to be stabilized to become a solvable problem, leading to the same expression for the bilinear form. There will be, however, a difference in right-hand-side terms. There is, however, no particular reason to choose one above the other since all formulations are consistent with perfect (i.e., divergence-free) measurements.

Remark 3 An alternative formulation can be introduced by re-defining Equation (3) as

∇ • w = -∇ • u meas + 1 |Ω| ∂Ω u meas • n, (6)
where the second term in the right-hand-side is needed to enforce the compatibility with respect to the boundary condition (4). However, as will be show later on, this leads w to have in general larger values than using the divergence-free model, even in places with no or little measurement errors, reducing its potential utility.

Remark 4

The regularity assumptions required to solve for w are stronger than the usual perturbations in real data, which are often not additive in nature and present strong jumps in space and time. Therefore, we do not aim that w is capable of "correcting" the data (e.g. by computing u meas + w) as it may be the purpose of e.g. data assimilation approaches. The field w rather aims to provide a complementary metric to detect defects in the measurements.

Remark 5

The choice of homogeneous Dirichlet boundary conditions is based on its consistency with the case of perfect measurements. In contrast, natural boundary conditions are not consistent: if w = 0 everywhere the pressure p cannot become the physical pressure on boundaries where homogenous Neumann boundary condition is enforced.

Remark 6 Note that the so-called "Stokes estimator" (STE) method for pressure reconstruction [START_REF] Švihlová | Determination of pressure data from velocity data with a view toward its application in cardiovascular mechanics. Part 1. Theoretical considerations[END_REF] is recovered dropping the first three terms of the left-hand-side.

Stabilized finite element formulation

The tetrahedral mesh with characteristic element size h obtained from the segmented medical image is denoted by Ω h , which is the discrete domain over we define the following functions spaces

V h = {w ∈ [H 1 0 (Ω h )] 3 : w ∈ [P 1 (K)] 3 ∀K ∈ Ω h } and Q h = {q ∈ L 2 0 (Ω h ) ∩ H 1 (Ω h ) : q ∈ P 1 (K)∀K ∈ Ω h }.
In order to aim for the clinical applicability, it is crucial to use fast and robust numerical schemes. For the spatial discretization, we adopt V h and Q h as spaces for w and p, respectively, using stabilized finite elements to ensure solvability. For the time discretization, we consider a backward Euler method with fixed time step τ to avoid GCL-type conditions. In order to avoid a root-finding problem at each time step the non-linear term on w will be treated semi-implicitly.

The resulting fully discrete stabilized formulation reads as follows. Given

w 0 = 0, for k ≥ 1 find (w k , p k ) ∈ V h × Q h such that B k (w k , p k ; v, q) = L k (v, q) (7)
for all (v, q) ∈ V h × Q h . The stabilized bilinear form is defined as:

B k (w, p; v, q) := A k (w, p; v, q) + S k conv (w; v) + S k press (w, p; v, q) (8)
with

A k (w, p; v, q) := Ω ρ τ w • v + ρ (u k meas + w k-1 ) • ∇ w • v + ρ w • ∇ u k meas • v (9) +µ∇w : ∇v -p∇ • v + q∇ • w
being the bilinear form associated to the non-stabilized weak form of (5), while the convection stabilization term is given by

S k conv (w; v) := Ω ρ 2 ∇ • (u k meas + w k-1 ) w • v
and the pressure stabilization term as

S k press (w, p; v, q) := δh 2 µ Ω h ρ (u k meas + w k-1 ) • ∇ w + ρ w • ∇ u k meas + ∇p • ρ (u k meas + w k-1 ) • ∇ v + ρ v • ∇ u k meas + ∇q
with δ > 0 some user-defined parameter. Finally, the right-hand-side is given by

L k (v, q) := Ω ρ τ w k-1 • v + k (v, q) with k (v, q) := Ω f k • v -µ∇u k meas : ∇v + δh 2 µ Ω h f k • ρ (u k meas + w k-1 ) • ∇ v + ρ v • ∇ u k meas + ∇q and f k = -ρ(u k meas -u k-1 meas )/τ -ρ u k meas • ∇ u k meas .
Additional theoretical properties of Problem [START_REF] Kozerke | Accelerating cardiac cine 3d imaging using k-t blast[END_REF] in terms of its solvability and time stability of the solution are shown in Appendix A.

Remark 7

The stabilization term S k conv is in general not consistent with the solution of (5). However, it is consistent when the measurements are perfect since, in such case w = 0 and ∇ • u meas = 0. The stabilization term S k press is weakly consistent due to the inclusion of h 2 .

3 Numerical experiments using synthetic data

Procedure for synthetic data generation

All testcases were created in the following way:

• A reference flow u is created using a finite element solver setup on an unstructured mesh of an aorta with a coarctation using high spatial and temporal resolutions. Physiologically relevant boundary conditions were set. The geometry with boundary conditions and resulting velocity field are shown in Figure 1. The simulation details are given in Appendix B. • The results of the simulation where interpolated in time with a time step of 0.03 s, hence downsampled by a factor of 30.

• In order to simulate the 4D Flow MRI acquisition, a complex magnetization field was computed for every component of the reference velocity as: M j = m 0 exp (iφ 0 + iπu j /venc) with j = 1, 2, 3, with u 1 , . . . , u 3 the time undersampled velocity fields. The value of φ 0 was assumed to be constant in space and time and equal to 7.5 • 10 -2 rad , and m 0 assumed constant and equal to 0.5 on the volume. The venc parameter setup will be explained in short. The background magnetization field was computed without any dependence on the velocity as M 0 = m 0 exp (iφ 0 ).

• Each magnetization field M 0 , . . . , M 3 was interpolated into a fine voxel-mesh with an element size of h b = 1 mm. Then, a smoothing filter on space via convolution with a Gaussian kernel with a standard deviation of 5h b was applied, in order to take into account partial volume effects. The filtered magnetization was interpolated into a grid mesh with the desired image resolution of 2 × 2 × 2 mm 3 around the reference mesh. Such interpolation was performed using a piecewise linear Lagrangian interpolator. Then, the four components were arranged as a multidimensional array, which will serve to create the images, denoted as M 0 , M 1 , M 2 , M 3 .

• We perturbed the magnetizations adding noise and aliasing by:

-The venc parameter was set lower than the peak true velocity. Therefore, when the velocity (in absolute value) exceeds the venc, the reconstructed velocity will be wrapped. Two venc parameters as the 150% and 70% of the maximum reference velocity were chosen, resulting in the values of 123 and 57 cm/s, respectively.

-The magnetization M j , j = 0, . . . , 3 includes Gaussian noise, having real and imaginary parts independent noise realizations with standard deviation 0.25. This results in a velocity noise with a variance of 17.73% of the maximum reference velocity, in the high venc case, and a variance of 11.41% in the lower venc case.

• Then, the 4D Flow measurements are given by

u meas = venc π angle( M j / M 0 ),
with / representing an element-wise division of the arrays and M j corresponding to the magnetization perturbed with the noise. Here, the time step index has been omitted in the description for the sake of readability.

• An image mask is created from the reference simulation on the uniform-grid mesh.

Then, a new semi-structured tetrahedral mesh following the aortic shape is created using the algorithm reported in [START_REF] Nolte | Validation of 4D flow based relative pressure maps in the aorta[END_REF]. The velocity is defined as a P 1 finite element field on such mesh to visualize the results and quantify the errors to the reference solution interpolated to the same mesh. The different velocity measurements generated are shown in Figure 2, first column from left to right.

Results

Problem (7) was solved using synthetic velocity measurements. In all the cases, δ = 1 was used for the pressure stabilization term, chosen as the smallest possible value regarding a local neighborhood insensitive to δ.

The results are shown in Figure 2 at the time of peak systole. For every data set, the two alternative formulations defined by Equations ( 3) and [START_REF] Baltes | Accelerating cine phase-contrast flow measurements using k-t blast and k-t sense[END_REF], refered now as Model A and Model B were tested, resulting in the fields w A and w B , respectively. Moreover, the difference field δu = u true -u meas is shown for comparison with w. Finally, the values for |∇ • u meas | are shown.

For the non-perturbed measurements (first row, (a)-(d)), w grows in zones where convective effects in the flow are more significant, i.e., mostly after the coarctation. It can be observed that when the measurements have no perturbations, w A is negligible compared with u meas . The exceptions are the convection-dominated regions where w achieves values below 30% of u meas . This could be attributed to the possibility that larger errors in the gradient of the measurements are more challenging to handle with the proposed stabilized finite element scheme. The divergence field shows similar behaviour. In contrast, w B presents very large values in a few spatial locations in spite that the measurements are perturbed only by the spatial blurring.

In the pure aliased case (second row, (e)-(h)), the divergence appears to be larger after the coarctation, where aliasing is present. However, in the ascending aorta it appears to detect aliasing only in a few places. This could be related by the fact that in the case of a "straight" flow, a wrapped velocity field does maintain the divergence. In contrast, w on both models considerably grows where aliasing occurs, compared to the non-perturbed case, but not as much as δu. Also, a coupling among the components of the field appears, therefore w on both models, does not point exactly as δu. Note that also not all points with aliasing are detected, most likely due to the homogenous Dirichlet boundary condition. In the case of w B , it turns out that its magnitude grows considerably in the regions where no aliasing is present, while the growth of w A does concentrate around the aliasing.

For the pure noise case, w B presents higher values than w A . Nevertheless, both present a similar "random" behaviour as δu but with lower magnitude in the case of w A due to its higher regularity. For the case with noise and aliasing, a combination of the two aforementioned behaviours is obtained. Note that since the velocity-to-noise ratio is proportional to the venc, in this case it is smaller than what is shown in the "pure noise scenario" with a larger venc. The divergence field shows again a small sensitivity to aliasing and exhibits an overall increment in the presence of noise.

Regarding the differences between both corrector fields, in summary w B presents larger values in artifact-free zones. w A , on the other hand, can better detect localized perturbations in the measurements. Moreover, as mentioned in Remark 5 and as it can be appreciated on Figure 2, the field w does not match δu due to its higher regularity than the perturbed measured velocities.

Consequently, since the purpose of w is to simply detect faults in the data, we recommend and adopt Model A for the rest of the manuscript.

Aditionally, a correlation study of the fields w and |∇ • u| for every case with respect to δu was performed. Only the results obtained by the Model A were used. Pearson and Spearman correlation coefficients were computed for each field at the moment of peak systole. Moreover, a rank of significant association was computed between the values of the fields using the Maximal Information Coefficient estimator (MIC e ) from [START_REF] Albanese | A practical tool for maximal information coefficient analysis[END_REF]. All these values are shown in Tables 1234. The w field shows better correlation with δu than |∇ • u meas | in cases where aliasing is present in the measurements. Without aliasing, for both fields the correlation is much lower, having w a slightly better performance. The MIC e coefficient shows higher level of significance between the fields when aliasing is present as well.

4 Numerical examples using real 4D Flow MRI data

Experimental phantom

A realistic thoracic aortic phantom was scanned using a clinical 1.5 T MR scanner (Philips Achieva, Best, The Netherlands) with a four-element phased-array body coil. The phantom was made of flexible silicone and a 11 mm orifice coarctation made of Technyl was placed in the descending aorta (for further details of the setup and the phantom see [START_REF] Urbina | Realistic aortic phantom to study hemodynamics using mri and cardiac catheterization in normal and aortic coarctation conditions[END_REF][START_REF] Montalba | Variability of 4d flow parameters when subjected to changes in mri acquisition parameters using a realistic thoracic aortic phantom[END_REF]). A blood mimicking fluid made with 60 % water and 40 % glycerol (Orica Chemicals, Watkins, CO) was used in the system. The fluid has a density of 1.119 g/cm 3 , dynamic viscosity of 0.0483 P and T1 value of 900 ms, which are in the range of values for human blood. The acquision was performed with a venc of 350 cm/s and using a cartesian sampling sequence with no k-space undersampling involved. In MRI, the noise level of the image increases when decreasing the voxel size. Therefore, three isotropic voxel sizes (coarse: 2.5mm, mid: 2.0mm and fine: 1.5mm) were acquired in order to investigate the results of w in front of different resolutions and SNR levels.

Figure 3 shows the 4D Flow measurements together with their result for w at the moment of peak systole. Also the divergence of the measurements are included.

First, note that for all three resolutions w tends to grow in zones where the flows becomes convection dominated, and not necessarily where the velocities are high, see Figure 3 (a) and (b). For instance, in the three regions of high velocity, the inflow tube and the descending aorta post stenosis, the inflow tube presents a negligible w value, where we expect The flow to be closer to a Stokes flow. This observation agrees with the results obtained in the synthetic case. However, real measurements have an additional source of error, namely the measurement technique assumes that the velocity field is constant within the time window of observation, in this case 40 ms. This induces to larger errors in the data that are not present in the synthetic case, possibly explaining the higher w values observed in the coarse voxel case.

When decreasing the voxel size, as expected from the reduction of the signal-to-noise ratio, w increases with the appearance of artifacts in the measurements. This occurs in the ascending aorta, in the stenosis, and in the inferior part of the arch. For all voxel sizes the divergence shows an increment post stenosis, with much lower values in the rest of the domain. Therefore, at least in this experiment the divergence of the measurements seems not to be capable of detecting spatial regions with decreased measurement quality, it only captures the regions where the velocity is larger.

Subjects

Two healthy volunteers were scanned in the previously described 1.5 Tesla scanner using a 4-channel torso coil. The local committee approved the study, and informed written consent was obtained from the participants. The acquisition parameters were: F OV in the range of 192 × 192 × 162 and 224 × 224 × 162 mm 3 , voxel resolution of 2.0 × 2.0 × 2.5 mm 3 , temporal resolution dt = 34 ms, venc = 150 cm/s, 25 cardiac phases, flip angle of 6 o and T R/T E = 4.9 ms/2.9 ms. From the resulting data, only the aorta was segmented in order to apply our in-house mesh generation algorithm. Afterwards, the problem was solved assuming a blood density of ρ = 1.2 gr/cm 3 and a dynamic viscosity of µ = 0.035 P , same values taken for the simulation with synthetic data. A study of the impact of these parameters on the solution w is presented in Appendix C.

Figure 4 shows the 4D Flow measurements, their divergence and the resulting w at peak systole. In all the cases w grows when flow artifacts appear in the measurements.

Volunteer 1 shows a highly regular blood flow which results in a small w. The divergence presents no large peaks. This suggests that the measurements are mostly perturbed by noise. Only in the ascending aorta w takes slightly higher values, probably due to larger convective effects in that region. The divergence field shows a concentrated spike on the top of the aortic arch likely to be caused by a boundary artifact, which is not detected by w due to the homogenous boundary condition.

Volunteer 2 shows a case where the velocity measurements are more perturbed, in particular on three different locations. In the ascending aorta, closer to the heart, the measurements presents discontinuities at certain locations. This is confirmed by the high values of w, while the divergence seems almost not to be perturbed, compared to other regions. In the aortic arch, all three fields appear to be perturbed. A last small region with larger w and divergence values appear in the distal part of the descending aorta.

Conclusion

We presented a new mathematical model -including a tailor-made discretization -to detect imperfections in full-field velocity measurements as 4D Flow MRI can obtain it. The derived Synthetic data experiments show the ability of the proposed approach to detect typical artifacts in 4D Flow MRI images, such as noise and aliasing, more robustly than computing the divergence of the measurements. With real 4D Flow MRI, we showed how w increases as the data quality decreases. That indicator mark regions where the data could be misrepresenting the blood flow, which is a valuable information when further flow quantification needs to be performed. In our experiments, we observed how some of these errors were not detected using the divergence as indicator.

At representative spatial and temporal resolutions of the data, computations of the proposed finite element discretization results take about 50 seconds without any parallelization using standard personal computers. That makes its clinical application feasible, since most of the burden lies in the blood lumen segmentation, where several flow markers are quantified.

Our work still presents some limitations. Concerning the numerical scheme, more advanced discretizations could be investigated to reduce the vector values arising from the discretization itself, particularly in convection-dominated regions. Moreover, the proposed discretization of w can not fully capture strong discontinuities in the measurements, e.g., in aliasing-contaminated measurements, mainly due to the regularity of the solution imposed by the discrete spaces used in this work. Therefore future work could consist in investigating the application of Discontinuous Galerkin approaches.

Another limitation of this study concerns the measurement generation in the synthetic data example. The analysis did not consider variability, due to a time blurring effect after several cardiac cycles, which occurs in the MRI measurement process due to sequential filling of the frequency space, respiratory motion artifacts, among others.

A Well posedness analysis of the discrete solution A.1 Existence and uniqueness of the discrete solution

The purpose is now to determine if problem ( 7) is well-posed by verifying if it fullfills Lax-Milgram theorem.

First, note that u k meas ∈ [H 1 (Ω h )] 3 for all k since the velocity data is bounded and interpolated to P 1 finite elements.

Denote the space of the whole solution vector W = V h × Q h , then we can prove the following results.

Lemma 1 The operator •

W : W → R defined by (v, q) 2 W := β v 2 H 1 (Ω h ) + S k press (v, q; v, q) is a norm on W = V h × Q h for β > 0.
Proof The first term defines a norm in V h and the second term is a seminorm in W , so • W is a seminorm in W . It remains to proof that (v, q) W = 0 ⇒ (v, q) = 0. Indeed, the first two terms have to be zero and therefore v = 0. Therefore, ∇q is a constant and also zero-valued since q ∈ Q h . Proposition 1 For v, w ∈ R d , C ∈ R d×s the following relation holds:

w • Cv ≤ C ∞ w 1 v 1 ≤ d C ∞ w 2 v 2 with m denoting the m -norm in R d .
Proof Directly from the proposition's definition.

Lemma 2 There exists α > 0 such that:

B k (w, p; w, p) > α (w, p) 2 W (10)
∀ (w, p) ∈ W \{0} under the condition:

ρ/τ + C -2 Ω µ/2 -ρ3 ∇u k meas ∞ > 0. (11)
Proof Using standard arguments, Poincaré's inequality and Lemma 1 the following relation holds:

A k (w, p; w, p) + S k conv (w; w) = Ω ρ τ w 2 2 + µ ∇w 2 2 + ρ(w • ∇)u k meas • w ≥ ρ τ + µ 2C 2 Ω w 2 L 2 (Ω h ) + ρ Ω (w • ∇)u k meas • w + µ 2 ∇w 2 L 2 (Ω h ) ≥ ρ τ + µ 2C 2 Ω -ρ3 ∇u k meas ∞ w 2 L 2 (Ω h ) + µ 2 ∇w 2 L 2 (Ω h ) ≥ β w 2 H 1 (Ω h ) with β = min(ρ/τ + C -2 Ω µ/2 -ρ3 ∇u k meas ∞ , µ 2
) > 0 under condition [START_REF] Mura | Enhancing the velocity data from 4d flow mr images by reducing its divergence[END_REF]. By adding the remaining term S k press (w, p; w, p), relation [START_REF] Thekla H Oechtering | Time-resolved 3-dimensional magnetic resonance phase contrast imaging (4d flow mri) reveals altered blood flow patterns in the ascending aorta of patients with valve-sparing aortic root replacement[END_REF] follows directly from the definition of the norm.

Lemma 3 There exists a constant M > 0 such that:

|B k (w, p; v, q)| ≤ M (w, p) W (v, q) W for all (w, p), (v, q) ∈ W . Proof Since |B k (w, p; v, q)| ≤ |A k (w, p; v, q)| + |S k conv (w; v)| + |S k press (w, p; v, q)|,
using Cauchy-Schwarz inequality and adding the missing norm terms in w we obtain

|S k press (w, p; v, q)| ≤ (w, p) W (v, q) W .
For the other terms, we can integrate the convective term by parts (hence canceling out the convective stabilization) and the pressure and divergence terms. Proceeding then using Cauchy-Schwarz:

|A k (w, p; v, q) + S k conv (w; v)| ≤ ρ τ + ρ ∇u k meas ∞ w L 2 (Ω h ) v L 2 (Ω h ) +ρ u k meas + w k-1 ∞ w L 2 (Ω h ) ∇v L 2 (Ω h ) +µ ∇w L 2 ∇v L 2 + ∇p L 2 (Ω h ) v L 2 (Ω h ) + ∇q L 2 (Ω h ) w L 2 (Ω h ) ≤ ρ τ β + ρ β ∇u k meas ∞ + ρ β u k meas + w k-1 ∞ + µ β (w, p) W (v, q) W + 1 √ β ∇p L 2 (Ω h ) (v, q) W + 1 √ β ∇q L 2 (Ω h ) (w, p) W .
And to end the proof we need to bound the pressure gradient:

∇p L 2 (Ω h ) ≤ ρ (u k meas + w k-1 ) • ∇ w + ρ w • ∇ u k meas + ∇p L 2 (Ω h ) + ρ (u k meas + w k-1 ) • ∇ w + ρ w • ∇ u k meas L 2 (Ω h ) ≤ µ δh 2 (w, p) W + ρ(u k meas + w k-1 ) L 2 (Ω h ) ∇w L 2 (Ω h ) + ρ∇u k meas L 2 (Ω h ) w L 2 (Ω h ) ≤ µ δh 2 + 1 √ β ρ(u k meas + w k-1 ) L 2 (Ω h ) + 1 √ β ρ∇u k meas L 2 (Ω h ) (w, p) W .
Lemma 4 There exists a constant C > 0 such that:

|L k (v, q)| ≤ C (v, q) W .
Proof We proceed using Cauchy-Schwarz inequality and adding the reminder terms to obtain the W -norm:

|L k (v, q)| ≤ ρ τ w k-1 L 2 (Ω h ) v L 2 (Ω h ) + | k (v, q)| ≤ ρ τ √ β w k-1 L 2 (Ω h ) (v, q) W + | k (v, q)|.
The last term can be bounded as:

| k (v, q)| ≤ f k L 2 (Ω h ) v L 2 (Ω h ) + µ ∇u k meas L 2 (Ω h ) v L 2 (Ω h ) + δ h 2 µ f k L 2 (Ω h ) δ h 2 µ ρ (u k meas + w k-1 ) • ∇ v + ρ v • ∇ u k meas + ∇q L 2 (Ω h ) ≤ 1 β f k L 2 (Ω h ) + µ β ∇u k meas L 2 (Ω h ) + δ h 2 µ f k L 2 (Ω h ) (v, q) W .
Theorem 1 There exists a unique solution in W of Problem [START_REF] Gamper | Compressed sensing in dynamic mri[END_REF] under condition [START_REF] Mura | Enhancing the velocity data from 4d flow mr images by reducing its divergence[END_REF] for all k > 0.

Proof Since w 0 = 0 ∈ [H 1 (Ω)] 3 , the bilinear and linear forms fulfill the requirements of the Lax-Milgram Theorem (i.e. Lemmas 2-4) for k ≥ 1.

A.2 Time stability of the discrete solution

We can furthermore prove the following energy balance:

Theorem 2 For (w k , p k ) solution of Problem (7), with (v, q) = 0 it holds w k 2 L 2 (Ω) ≤ w k-1 2 L 2 (Ω) (12)
under the condition

µ ≥ C 2 Ω 3ρ ∇u k meas ∞ (13)
Proof Testing [START_REF] Gamper | Compressed sensing in dynamic mri[END_REF] with v = w k and q = p k and using similar arguments as in Lemma 2, it is obtained

ρ 2τ w k 2 L 2 (Ω h ) - ρ 2τ w k-1 2 L 2 (Ω h ) = -µ ∇w k 2 L 2 (Ω h ) - Ω ρ w k • ∇ u k meas • w k - ρ 2τ w k -w k-1 2 L 2 (Ω h ) -S k press (w k , p k ; w k , p k ),
where the two first terms in the right-hand-side come from the continuous problem, and the two last terms are dissipative due to the numerical scheme. Bounding the former using Poincaré's inequality and the later by zero we obtain

ρ 2τ w k 2 L 2 (Ω h ) - ρ 2τ w k-1 2 L 2 (Ω h ) ≤ (-µC -2 Ω + 3ρ ∇u k meas ∞ ) w k 2 L 2 (Ω h ) ,
which combined with condition (13) leads to relation [START_REF] Busch | Reconstruction of divergence-free velocity fields from cine 3d phase-contrast flow measurements[END_REF].

Remark 8 Condition (13) is very conservative and in the test cases we observe that the condition is not satisfied. However, in our computations, we have never obtained instabilities, even in situations with temporal jumps in the data as when velocity aliasing appears.

B Reference synthetic aortic flow

Synthetic data was generated in a drawn-aortic mesh with a mild coarctation in the descending aorta and consist of 2,752,064 tetrahedrons and 510,755 vertices. An incompressible Navier-Stokes problem was solved with a monolithic formulation with P 1 /P 1 finite elements and with Temam and PSPG stabilizations. The boundary conditions were set as follows.

At Γ inlet the following velocity profile was imposed:

u inlet = u stokes f (t) n,
where u stokes corresponds to a profile resulting from solving a transient Stokes problem with a non-homogenoeus Neumann boundary condition on Γ inlet in the same geometry, while f (t) is a given waveform defined as

f (t) =    sin( πt T ) if t ≤ T π T (t -T )e -κ(t-T ) if T c > t > T.
Here, T gives the opening time of the heart valve, T c the total duration of the cardiac cycle and 1/κ represents the typical time for the closing of the valve. Then, a non-slip condition for the velocity at the walls Γ wall was imposed. Finally, a three-element Windkessel boundary condition for the rest of the outlets was applied, i.e.: µ ∂u ∂n -pn = -P (t)n on Γ w, , = 1, . . . , K

where K denotes the number of outlets. The term P (t) is computed by solving the equations:

             P = R p, Q + π Q = Γ u • n C d, dπ dt + π R d, = Q , (14)
where R p, and R d, represent the resistance of the vasculature proximal and distal to Γ , respectively, and C d, the compliance of the distal vessels. The Windkessel constants (R p , C, R d ) were tuned by hand in order to have a standard physiological flow regime to achieve approximately 70%/30% split in the peak flow rate between the descending aorta and supra-aortic branches [START_REF] Urbina | Realistic aortic phantom to study hemodynamics using mri and cardiac catheterization in normal and aortic coarctation conditions[END_REF][START_REF] Alastruey | On the impact of modelling assumptions in multi-scale, subject-specific models of aortic haemodynamics[END_REF].

For the physical parameters of the fluid and the constants of the Windkessel models see Table 5. Additionally, blackflow stabilization was added in every outlet of the system [START_REF] Bazilevs | Patient-specific isogeometric fluid-structure interaction analysis of thoracic aortic blood flow due to implantation of the jarvik 2000 left ventricular assist device[END_REF]. The initial conditions are set as u 0 = 0 and π 0 = 85 mmHg for = 1, . . . , K, which correspond to approximately the periodic state of the 3D-0D system. The simulation was performed with a total time of T f = 0.8 s with a timestep of dt = 0.001 s.

Parameter

C Viscosity sentivity study

In order to know how sensitive is the corrector field w to variations of the dynamic viscosity parameter µ, we computed two fields using the values µ -= 0.0175 P and µ + = 0.070 P with the 4D Flow data set acquired for the healthy volunteers. These values are in a 50% range of reference value for human blood used in this work (µ 0 = 0.035 P ), more concretely, we computed the H 1 -norm time curves of the corrector field over time, which are shown in Figure 5. Note that, as the second volunteer has relatively higher velocities, the effect of the viscosity will be "hidden" by the dominance of the convective forces within the fluid. It can be seen that higher differences start to appear during diastole, when the Reynolds number are relatively small compared with the starting of the cycle. This can be explained by the dominance of drag forces in a low velocity regime. In Volunteer 1, the relative difference of the H 1 -norm is no more than a 15% of the values obtained using µ 0 . The latter can be confirmed from the Figure 6, in which the fields are shown at t = 0.6 s, corresponding to the time instant where differences are maximal. In conclusion, within a reasonable interval of viscosity values no major differences in the corrector fields are observed. Therefore, we can safely adopt the reference value µ 0 = 0.035 P .

(a) w using µ -for volunteer 1 (b) w using µ 0 for volunteer 1 (c) w using µ + for volunteer 1 (d) w using µ -for volunteer 2 (e) w using µ 0 for volunteer 2 (f) w using µ + for volunteer 2

Figure 1 :

 1 Figure 1: Boundary conditions for the simulation with the reference solution obtained for generate synthetic data sets. In (a) boundary conditions, (b) Velocity vector field in the reference mesh at peak systole.

Figure 2 :

 2 Figure 2: Measurement, w A,B , error and divergence field for different measurements: (a)-(e) non perturbed , (f)-(j) aliased, (k)-(o) noisy and (p)-(t) aliased and noisy.

Figure 3 :

 3 Figure 3: Results for aortic phantom data, for three spatial resolutions of the measurements: measurements (left), (mid )w, and (right) divergence of the measurements.

(a) u meas for volunteer 1 (Figure 4 :

 14 Figure 4: 4D Flow measurements, the corresponding w, velocity and divergence fields for the volunteers at the time of peak systole.

4 Table 5 :

 45 Value ρ (gr • cm 3 ) 1.2 µ (P ) 0.035 U (cm • s -1 ) 75 T c p (dyn • s • cm -5 ) 480 520 520 200 R d (dyn • s • cm -5 ) 7200 11520 11520 4800 C (dyn -1 • cm 5 ) 4 • 10 -4 3 • 10 -4 3 • 10 -4 4 • 10 -Physical parameters and numerical values of the three-element Windkessel parameters for every outlet.

  (a) ||w|| H1 for volunteer 1 (b) ||w|| H1 for volunteer 2

Figure 5 :

 5 Figure 5: Comparison between norms of the corrector field under three different dynamic viscosity values.

Figure 6 :

 6 Figure 6: Corrector fields at t = 0.6 s. The arrow scale have been remained the same with respect to the scaling used in the rest of this work for a better comparison.

Table 4 :

 4 Noisy + aliased data

		w	|∇ • u meas |		w	|∇ • u meas |
	Pearson coeff. 0.1493	-0.0186	Pearson coeff. 0.5939	0.2023
	Spearman coeff. 0.2214	0.0588	Spearman coeff. 0.2510	0.1258
	MIC e	0.0366	0.0141	MIC e	0.2540	0.1292
	Table 1: Non perturbed data	Table 2: Only aliased data
		w	|∇ • u meas |		w	|∇ • u meas |
	Pearson coeff. 0.0026	0.0035	Pearson coeff. 0.4911	0.2045
	Spearman coeff. 0.0080	-0.0002	Spearman coeff. 0.1534	0.0681
	MIC e	0.0060	0.0053	MIC e	0.1303	0.0360
	Table 3: Only noisy data