Geometrically induced wall shear stress variability in CFD-MRI coupled simulations of blood flow in the thoracic aortas

Aortic aneurysm is associated with aberrant blood flow and wall shear stress (WSS). This can be studied by coupling magnetic resonance imaging (MRI) with computational fluid dynamics (CFD). For patient-specific simulations, extra attention should be given to the variation in segmentation of the MRI data-set and its effect on WSS. We performed CFD simulations of blood flow in the aorta for ten different volunteers and provided corre- sponding WSS distributions. The aorta of each volunteer was segmented four times. The same inlet and outlet boundary conditions were applied for all segmentation variations of each volunteer. Steady-state CFD simula- tions were performed with inlet flow based on phase-contrast MRI during peak systole. We show that the commonly used comparison of mean and maximal values of WSS, based on CFD in the different segments of the thoracic aorta, yields good to excellent correlation (0.78 – 0.95) for rescan and moderate to excellent correlation (0.64 – 1.00) for intra- and interobserver reproducibility. However, the effect of geometrical variations is higher for the voxel-to-voxel comparison of WSS. With this analysis method, the correlation for different segments of the whole aorta is poor to moderate (0.43 – 0.66) for rescan and poor to good (0.48 – 0.73) for intra- and interobserver reproducibility. Therefore, we advise being critical about the CFD results based on the MRI segmentations to avoid possible misinterpretation. While the global values of WSS are similar for different modalities, the variation of results is high when considering the local distributions. the 4D-Flow MRI segmentations. we show that the variability is for both rescan and coupling shows the potential for studying the progression of in serial follow up scans.


Introduction
Aorta pathologies such as an aortic aneurysm, dissection, and coarctation affect its function [1]. These conditions are not completely understood. However, aberrant blood flow patterns and associated hemodynamical parameters like wall shear stress (WSS) have been linked to aortic dysfunction [2,3].
To study WSS in the aorta, we need to have information about the blood flow and the location of the aortic wall. In theory, both of these can be acquired using 4D-flow magnetic resonance imaging (MRI) [4]. While the 4D-flow MRI acquisition can produce a velocity field that closely represents the actual aortic flow [5][6][7], using this technique to acquire the location of the aortic wall may be incorect [8,9]. The irregularities in the segmented wall location can be influenced by the movement of the aorta during the cardiac cycle as well as due to its compliance [9,10]. Besides, the segmentation of the geometry is done (partly) manually in most commercial and research tools. Therefore, the observer's lumen interpretation results in segmentation variability [10]. Van der Palen et al. [9] showed that this variability, among others, may lead to differences in MRI estimated WSS as high as 30%. Additionally, while several methods have been proposed to calculate WSS from the 4D-flow MRI velocity field [11][12][13][14], the MRI-based WSS is underestimated due to the resolution of the flow field [3,15]. Hence, to be able to accurately predict WSS, a different method should be used.
To do this, many studies showed that coupling MRI with computational fluid dynamics (CFD) leads to patient-specific results that closely represent reality [16][17][18]. An important aspect, that is often not addressed, of evaluating the accuracy of patient-specific simulations is the influence of variability in geometry reconstruction from imaging techniques on numerical results. The effects of geometry variation on the calculated blood flow and associated parameters are assessed in two CFD studies focused on coronary vessels [19,20]. Both concluded that a small variation in geometry has a considerable effect on the predicted hemodynamics. However, these conclusions were based on just one case and hence the wide applicability is questionable. Moreover, none of these studies focused on the aorta, where the biggest variation in geometry can be expected due to its extensive motion during the cardiac cycle.
The present study aims to bridge this gap in literature and evaluate the effect of variations in the input geometry on WSS obtained from CFD. By this, we extend the previously performed 4D flow MRI based segmentation variability [10] and WSS variability [9] and we aim to establish a baseline for variability in CFD simulations. We perform simulations on aortas from ten volunteers at the peak systole. Four different segmented geometries are available for each volunteer. We compare the simulated WSS using statistical analysis as shown in Ref. [9], where the effect of varying geometry on WSS is evaluated based on the comparison of average values in five investigated regions -zero-dimensional analysis. Moreover, we extend this approach and propose two different ways of higher-dimensional comparison. First, we compare WSS between the different geometries based on its circumferential average alongside the centerline -one-dimensional analysis. Second, we a perform voxel-to-voxel comparison of flattened surface maps -two-dimensional analysis. We show the influence of the method for analysis on the agreement, where the WSS correlation between the different geometries decreases with the increasing order of the method.

Methodology
The present study is based on the data presented in our previous studies [9,10], where ten healthy volunteers were analyzed, but is now extended with CFD simulations. The characteristics of the volunteers, like gender, age, weight, height, and heart rate, can be found in Ref. [9] and the morphometric characteristics of the acquired aortas can be found in Ref. [10]. All volunteers gave informed consent before the MRI acquisition and METC approval was obtained for the study of these volunteers.

MRI data-set
The studied data set consists of ten different volunteers with four different geometries for each. 4D-flow MRI was performed on a 3.0 T Scanner (Ingenia, Philips Medical Systems, Best, The Netherlands with Software Stream 4.1.3.0). The spatial resolution was 2.5 × 2.5 × 2.5 mm 3 and the temporal resolution was 35.1-36.5 ms. Additional details about the aortic 4D acquisition can be found in Ref. [9]. All volunteers were scanned twice, with a 30-min interval between the scans.
Afterward, acquired 4D-flow MRI data were segmented using CAAS MR 4D flow v1.1 (Pie Medical Imaging BV, Maastricht, The Netherlands). A semi-automated segmentation algorithm was used that optimizes the location of the surface while maintaining the smoothness of the 3D surface [9]. The segmentations were manually adapted from the level of the aortic valve till the thoracic descending aorta (at the same level of the aortic valve).
The four different geometries of each volunteer were obtained with different segmentations, as follows: (i) 4D-Flow MRI was performed and the geometry was segmented by the first observer -RvdP (geometry Ascan), (ii) the MRI scan was reproduced with 30 min from the first scan with consequent segmentation by RvdP (geometry B -rescan), (iii) the segmentation by RvdP was reproduced on the data from the first scan (geometry C -intraobserver), (iv) the last geometry was segmented from the first MRI data set by the second observer -PB (geometry D -interobserver). The repeated analysis was performed blindly to the results of the previous analyses. RvdP has six and PB has fifteen years of experience in cardiovascular MRI. All geometries are shown in Fig. 1. In this study, only the thoracic aorta without the branches is evaluated. WSSMRI was estimated from the MRI velocity field using CAAS MR 4D flow v1.1 (Pie Medical Imaging BV, Maastricht, The Netherlands). For each wall-point of the segmented surfaces, WSS was calculated using a quadratic approximation of the axial velocity profile perpendicular to the aortic wall.
For analysis purposes, the resulting surfaces were split into five different investigated segments using cut-planes based on anatomical landmarks, using CAAS MR 4D flow v1.1: sinotubular junction (1), midascending aorta (2), origin of the innominate artery (3), beyond the left subclavian artery (4), the mid-descending thoracic aorta (5), and the descending aorta at the level of the aortic valve (6). This resulted in five segments (from the sinotubular junction to the descending aorta): proximal ascending aorta (pAAo), distal ascending aorta (dAAo), aortic arch (Arch), proximal descending aorta (pDAo), and distal descending aorta (dDAo). The segments with the corresponding cut planes are shown in Fig. 2a.

Geometry preparation
The raw geometry obtained directly from MRI is not suitable for performing CFD simulations, and hence it has to be pre-processed. Several improvements of the surface have to be performed such as smoothing, capping ends, and the addition of flow extensions at the outlets and inlets. The exact parameters for the surface manipulation using Vascular Modeling Toolkint (VMTK) [21] are as follows: • smoothing -Taubin method, passband 0.4 • surface subdivision -butterfly method • flow extensions -adaptive length of 3R i , where R i is the mean profile radius of the corresponding inlet and/or outlet The value for smoothing was chosen to obtain a smooth surface while maintaining its size. An example of the original MRI data set and processed geometry is shown in Fig. 2.

Discretization
An example of a computational mesh is shown in Fig. 2b. The meshing is done using ICEM CFD (Ansys Inc., Canonsburg, PA, United States). For every geometry, a hybrid numerical mesh was constructed with a fine boundary layer. The boundary layers are covered by prismatic control volumes, whereas the tetrahedral elements are applied in the center region of the aorta. Note that this refined numerical mesh in the proximity of the vessel wall is applied to be able to properly capture velocity gradients within relatively thin boundary layers. An order of magnitude estimate for usual flow in aorta yields the boundary layer to be δ ∼ 0.1 mm. Therefore, the cells close to the wall should be at least of this size [22].
The boundary layer was constructed using ten layers, an exponential growth rate of 1.2, and constant first-layer thickness of 0.01 mm. The computational mesh used for the simulations contained between three and five million control volumes. This resolution proved to be sufficient to obtain the grid-independent solutions (details can be found in appendix B).
Besides, Fig. 2 shows the cut planes which were used to divide the geometry for the analysis. Only segments between the cut-planes were used for the analysis. This division is important since the original geometries (shown in Fig. 1) vary in the arterial length due to different segmentation.

Governing equations and simulation setup
The steady blood flow is described by the conservation of mass and momentum (Navier-Stokes) equations for the incompressible fluid, which are written as: where u is the velocity vector, ρ is the density (ρ = 1060 kg/m 3 ), p is the pressure and μ is the dynamic viscosity (μ = 3.5 mPa⋅ s). Flow was assumed to be laminar, with peak Reynolds number varying from 2024 to 3142. The boundary condition for the inlet face was set to a mass flow inlet with a flat profile and the value based on the MRI measurements. The mass flow rate was exported from the reformatted 4D flow measurements and was kept the same for all four segmentations per volunteer. The outlet was set to the outflow boundary condition, with a zero diffusion flux for all flow variables. The rigid wall assumption with the no-slip velocity boundary condition was imposed at the aorta wall. All simulations were performed using Ansys Fluent 18.2 (Ansys Inc., Canonsburg, Pennsylvania, USA). The details of the numerical solver setup are listed below: • Solver -pressure based • Pressure-Velocity Coupling -SIMPLE • Spatial discretization  -Gradient -Least Squares Cell-Based -Pressure -Second-Order -Momentum -Second-Order Upwind • Residuals -10 − 5

Post-processing
To evaluate the variation between the calculated WSS in the different geometries obtained from MRI, a two-dimensional mapping approach is performed for all three-dimensional aortic walls. In this approach, the contours of WSS at the aorta wall are divided into several circumferential and longitudinal sections. In total, 100 circumferential segments were created for each surface and the length of each segment in the longitudinal direction was 0.5 mm. The surface is then cut along the inner curvature of the aorta wall, making a simple two-dimensional projection map of the WSS distribution. The most important steps in this mapping approach are shown in Fig. 3. Note that for illustration purposes, in this particular case, segments were significantly bigger (10 circumferential sectors with 1 cm in the longitudinal direction). In the last step, the 2D maps of the WSS from different segmentations were superimposed using an in-house image registration protocol.

Statistical analysis
Statistical analysis was performed on the acquired data using an inhouse code implemented in Matlab R2018b (MathWorks, Natick, MA, USA). The WSS mean and WSS max are presented together with the standard deviation. The acquired data were analyzed in three different ways.
First, the aorta wall surface was divided into five sections and the mean and maximal value of WSS were calculated for each section, resulting in a single WSS parameter (consequently, we named this approach -'Point Analysis'). The Bland-Altman analysis was performed on these data, where the mean difference and the limits of agreement were calculated in the mean and max WSS, respectively, between the different surfaces. The limits of agreement are equal to 1.96σ (where σ is the standard deviation). Moreover, the Spearman Correlation Coefficient (ρ) was calculated for both the mean and maximal values of WSS. This analysis was equal to the one based on the MRI data of [9].
Secondly, the circumferential average and circumferential maximal values were analyzed alongside the whole centerline of the studied aortas (we name this approach -'Line Analysis'). These data were subsequently analyzed similarly as the 'Point Analysis'. Plots of the mean difference in the circumferential average or circumferential maxima were created together with the limits of agreement. Also, ρ calculated for the five different regions of the thoracic aorta.
Lastly, the statistical analysis was performed on the created WSS maps (we name this approach -'Surface Analysis'). For the surfaces, WSS mean and the standard deviation from the different surfaces were calculated for each of the volunteers. Afterward, the voxel-to-voxel ρ was calculated for all volunteers based on the whole surface, as well as the five investigated regions of the thoracic aorta.
For all methods, the classification of ρ was as follows: (ρ ≥ 0.95: excellent); (ρ ∈ (0.95; 0.85〉: strong); (ρ ∈ (0.85; 0.70〉: good); (ρ ∈ (0.70; 0.50〉: moderate); (ρ < 0.5: poor). The significance level was set to p = 0.05, hence the analysis with p-value p < 0.05 is considered as statistically significant. Tables 1 and 2 show the results from the Bland-Altman analysis [23] for the mean and maximal values of WSS, respectively, of the investigated regions of the thoracic aorta. The analysis consists of the mean difference between the WSS for the different segmentations and corresponding limits of agreement.

Point analysis
The analysis based on the average WSS in the thoracic aorta shows a mean difference between the values in a range from 0 Pa to 1 Pa (if absolute values are taken into consideration). The limits of agreement lie between 1 Pa and 3.7 Pa. The results for WSS max in the thoracic aorta sections show higher values for both the mean difference and the limits of agreement. For this case, the mean difference lies between 0 Pa and 9.74 Pa (in absolute values) and the limits of agreement are from 6 Pa to 23 Pa.
The analysis of the MRI results performed by one observer (RvdP) [9] considers just the following surface combinations: A-B (scan-rescan), A-C (intraobserver), and A-D (interobserver). The results presented here show a similar trend to the ones presented by RvdP [9] The variability of WSS and the limits of agreement are higher for regional WSS max . However, the results obtained from CFD show higher mean differences Next, the Spearman Correlation Coefficient was calculated for both the mean and maximal values of WSS in the investigated regions of the thoracic aorta. The results are presented in Table 3. For WSS mean the correlation between the surfaces is good to strong for most of the cases. The only exceptions are surface combinations A-D in the pAAo where the correlation is moderate and A-D in the pDAo where the correlation is excellent. When looking at the average values per region, the correlation increases from the pAAo (good) reaching its maximum in the pDAo (strong) and slightly decreasing in the dDAo (good). Not many differences can be found between the different surface combinations.
The maximal values of WSS show more variety, with the correlation being moderate to strong. However, a case with excellent correlation can be found (C-D in pDAo) as well as cases with poor correlation (A-C and B-C both in dDAo). Unlike for the WSS mean , the correlation does not reach its maximum in the pDAo but rather gradually decreases. The average of respective correlation for all surface representations decreases from good in the pAAo to moderate in the dDAo. This trend can be explained by the overall values for WSS. From Table 2 can be seen, that both the difference between the value of WSS max as well as the limits of agreement increase downstream the aorta. A similar trend can be seen also in the actual values of WSS for each of the studied surfaces. WSS is gradually increasing from the aortic root.
From the paper by Van der Palen et al. [9], the results show the correlation to be good to strong for WSS mean and moderate to good for WSS max in the case of A-B. For the surface combination A-C the correlation was found to be good to excellent for both WSS mean and WSS max . Finally, results for A-D show good to excellent correlation for both WSS mean and WSS max .
The results from the CFD analysis show a similar correlation of mean WSS for all cases and a lower correlation of maximal values of WSS. The difference in WSS max arises especially in the descending part of the thoracic aorta. While the correlation based on MRI was overall good to excellent [9], the agreement for CFD is moderate to good in these regions.
An example of the plots for the scan-rescan analysis (A-B) of WSS mean and WSS max based on the Bland-Altman analysis is shown in Fig. 4 together with the correlation plots. Plots for the other surface comparisons can be found in appendix D.

Line Analysis
In the next step, the circumferential averages of WSS were calculated along the corresponding centerline for all geometries. Afterward, the Table 1 Bland-Altman analysis for WSS mean showing mean difference and the limits of agreement (1.96 times standard deviation σ) in the investigated regions of the thoracic aorta: proximal ascending aorta (pAAo), distal ascending aorta (dAAo), aortic arch (Arch), proximal descending aorta (pDAo), and distal descending aorta (dDAo).  Table 2 Bland-Altman analysis for WSS max showing mean difference and the limits of agreement (1.96 times standard deviation σ) in the investigated regions of the thoracic aorta: proximal ascending aorta (pAAo), distal ascending aorta (dAAo), aortic arch (Arch), proximal descending aorta (pDAo), and distal descending aorta (dDAo).  average of WSS mean per volunteer was calculated together with the standard deviation. The resulting plots for all ten volunteers can be seen in Fig. 5. As shown in Fig. 5, the agreement between the geometries is better in the pAAo and dAAo than in the rest of the aorta. In the arch, the variation of WSS mean shows higher discrepancies as the flow starts to be more complex. Further downstream the aorta, in the pDAo and dDAo the agreement between WSS mean obtained from different segmentations, is very weak.
The Spearman correlation coefficient was calculated for the circumferential averages of WSS alongside the centerline. The values for all possible surface combinations are shown in Table 4. The correlation between the different surfaces is better for the circumferential averages. Mostly moderate and good correlation can be reported. The average among all volunteers shows a good correlation for A-C, A-D, and C-D and moderate correlation for A-B, B-C, and B-D. However, the differences are relatively small. The standard deviation varies between 0.05 and 0.11 and the p-value was smaller than 0.0001 for all cases.

Surface analysis
The results of the two-dimensional mapping procedure of WSS along the thoracic aorta wall are shown in Figs. 6 and 7. It can be seen that WSS is lower in the ascending part and higher in the arch and the descending aorta. In the descending aorta, significant variations are present in local WSS for all geometries. These results agree with the circumferential averages, where the variation between the WSS was prominent in the descending aorta and suggests that the possible error in WSS based on CFD will be most prevalent in these parts.
Until this point, the WSS distributions were assessed qualitatively by visual comparison. To evaluate more quantitatively the degree of variation in WSS caused by the geometry diversity, the map of average WSS and the standard deviation were additionally created from the four segmentations (shown as the last two inserts in Figs. 6 and 7). Two local maxima located between the aortic arch and the proximal descending aorta can be identified for all presented geometries. The plot of standard deviation shows the highest values in the descending part of the aorta as well.
The results of the Spearman correlation coefficient of WSS based on the voxel-to-voxel analysis of the whole surface representing the thoracic aorta can be found in Table 5. The results show that the voxelto-voxel correlation is mostly moderate, and in a few cases, good or poor. The average of all volunteers shows a moderate correlation for all possible surface comparisons. The standard deviation for the Spearman Correlation Coefficient is between 0.06 and 0.09 among the different segmentation combinations. It should be also noted, the p-value for all cases was smaller than 0.0001.
The correlation coefficients do not show a better or worse agreement among different combinations of surfaces. All segmentations made from the scan (A, C, D) have a similar agreement. Also, the segmentation based on the re-scan (B), which was performed in slightly different physiological conditions, has a similar agreement with the scan-based surfaces as can be seen from Table 5. Hence, the CFD based WSS is not highly affected by the re-scan conditions. The Spearman correlation coefficient based on the whole surface has lower values than the one based on the circumferential averages showed in Table 4. This can suggest that while the locality of WSS max and WSS min may differ, the mean values of the cross-sections are close to each other.
Finally, we have calculated the voxel-to-voxel Spearman correlation coefficient based on the investigated regions of the thoracic aorta (as was shown in Fig. 2). The results for the correlation coefficient in the distinct regions can be seen in Table 6. A similar trend can be observed for all surfaces, where the correlation gradually decreases from the inlet to the outlet of the aorta. The highest correlation can be found at the beginning of the aorta -in the pAAo. In this region, a moderate correlation can be found for all surface combinations. In the dAAo, the agreement is lower but still moderate. In the aortic arch, the correlation is even lower, on average still being moderate but with values around Table 4 Spearman correlation coefficient r based on the circumferential averages for ten volunteers (V1-V10) with the mean value and standard deviation (SD  These results confirm the findings in the comparison of the surfaces from Figs. 6 and 7 as well as the graphs of circumferential averages shown in Fig. 5. The agreement between the surfaces is lower in the pDAo and dDAo. It has to be noted that the agreement based on the voxel-to-voxel comparison is significantly lower in comparison to both the results observed for MRI by van der Palen [9] and the 'Point Analysis' for CFD. Fig. 7. Maps of WSS for volunteer 6-10 with the mean WSS map per volunteer and the corresponding standard deviation (SD) map; all of the maps are split in the five investigated regions: pAAo -proximal ascending aorta, dAAo -distal ascending aorta, Arch -aortic arch, pAAo -proximal descending aorta, pAAo -proximal descending aorta.

Table 5
Spearman correlation coefficient r based on the voxel-to-voxel analysis for ten volunteers (V1-V10) with the mean value and standard deviation (SD).

Discussion
Image-based CFD of arterial geometries is gaining popularity to study various arterial diseases. While numerous advances in simulating the flow in patient-specific arteries have been developed, many questions have also been raised about the validity of such simulations. There are many factors in the aortic flow that can have a significant impact on the generated flow patterns. These include the types of the inlet and outlet boundary conditions, movement of the wall, pulsating character of the inflow, blood rheology, etc. A great majority of these factors have been previously investigated, mostly on small cohorts.
However, the base of the simulations -the geometry of the aortaaffects the flow as well and the effects of small variations in geometry on the WSS, have often been overlooked in the MRI-based CFD studies. Particularly for MRI, it was previously shown that the morphometric characteristics of segmented aorta varies between rescan, intra-and interobserver analysis [10] and this has an effect on WSS estimated based on the 4D-flow MRI velocity fields [9]. Hence, it is of great importance to understand, how small variations in geometry due to the segmentation process affect the simulations. In this study, we aimed to answer this question and create a confidence interval for both the expected error as well as the method, with which such variability should be investigated.

Evaluation method
In MRI based studies, the WSS variability is often evaluated based on a commonly used 'Point Analysis', in which WSS mean or WSS max of a few segments of the aorta are compared between the different segmentations. This leads to many studies reporting a high level of agreement in WSS between intra/interobserver and scan/rescan [7,9,24].
We have performed the same point analysis also on our CFD data-set, which resulted in a good agreement between different surfaces, similar to MRI. While the agreement in the ascending aorta and arch are close to MRI, more differences are observed in the ascending aorta. The discrepancies between CFD and MRI variability based on the 'Point Analysis' are likely caused by higher velocities in this region due to the gradual narrowing of the artery and the omission of the branching arteries. Higher velocity causes higher WSS and steeper gradients close to the wall. Consequently, small irregularities between different segmentations eventually lead to a worse correlation between them. Additionally, for MRI, the velocity field does not change with different segmentations (except for scan/rescan comparison). For CFD, the velocity field may slightly differ due to the differences in the geometry, which also adds to the higher disagreement in the descending aorta.
While this type of analysis is widely accepted by the medical society, considering just a few data-points per segmentation and high spatialaveraging of the results may lead to over-estimation of the quality of agreement. By visual inspection of MRI-based WSS contours (Fig. A8b), multiple discrepancies between different segmentations can be observed, which are not represented in the results based on the 'Point Analysis'. This shows, that performing simplified analysis and comparing the results just in terms of WSS mean or WSS max in the particular segments of the aorta is insufficient since a lot of the information is lost due to averaging.
We observed an identical trend also in WSS CFD . While the 'Point Analysis' shows a high agreement, 'Line Analysis' and 'Surface Analysis' give on average moderate and poor agreement respectively. Both of these methods give much more information about the actual WSS trends in the aorta and can capture the local variations, which may be vital for evaluating disease progression. Additionally, the WSS correlation is lower in the descending part of the aorta as shown in the more detailed analyses. This was not visible from the analysis of MRI results nor from the 'Point Analysis' of the CFD results.

Effect of morphology on WSS
Morphology of aorta based on 4D flow MRI has a certain degree of uncertainty [10]. This also brings variability to CFD-based WSS. The three method of analysis show that the WSS variability gradually increases from the ascending aorta to the descending aorta. The higher value of WSS in the descending aorta is due to the smaller radius of the descending aorta compared to the ascending part [10]. Also, the flow features in this region are more complex. The complex flow features are originating from interactions of recirculation zone and secondary flow structures (so-called Dean vortices [25]). While all surfaces exhibit similar distributions, some differences can be observed as well, especially comparing the segmentation C to the rest. The maps of WSS from segmentations C were smoother than the others. This is reflected in the smaller variations in maximal/minimal WSS in all parts of the aorta compared to the other segmentations.
Accurate 4D flow MRI-based segmentations are more difficult to obtain since the method is highly dependent on the spatial and temporal resolution and the velocity field. This often leads to the method not being completely automatic and many manual adjustments need to be applied during the segmentation procedure. Consequently, the segmentations may vary, even for the same MRI data-set.
This affects the simulations and it is crucial that the error in WSS due to the segmentation variation is taken into account. Whereas the effects of the different segmentations are not high and the global WSS mean and WSS max correspond well, many differences can be seen in the details of WSS. Similar results were presented also for intracranial aneurysm [26, ], where WSS varied between 28% and 51%, depending on segmentation. Because of this, a high-quality segmentation is desirable for an accurate estimation of WSS.
Additionally, as was shown for intracranial aneurysm, the branch positioning and diameter have a great effect on WSS [26,]. It should be noted that this study omits the main branching arteries in the aortic arch. Consequently, the flow rate in the descending aorta is higher for CFD than for the MRI. Including the branches in the segmentation would probably lead to even higher discrepancies, since both the positioning of the branching arteries as well as their diameter could be different for four segmentations. This, in combination with more complicated boundary conditions for the branching arteries, could result in more variation in WSS between segmentations. This stresses the importance of high-quality segmentation for accurate estimations of WSS.

Clinical implications
WSS is a patho-physiological stimulus at the intimal surface of the aortic vessel wall that has been shown to alter gene expression and endothelial cell function [27]. Altered shear stress, either in longitudinal or circumferential direction, can promote endothelial changes that can create an area at risk for vascular remodeling and aneurysm growth. Accurate mapping of the variations in wall shear stress may prove to become a very relevant clinical tool, but before its introduction, the effect of accurate aortic lumen segmentation and observer-and repeated scanning-induced variations on the wall shear stress quantitation needs in-depth evaluation. In this study we have used combined MRI-CFD approach to study the variations on wall shear stress.
Our analyses showed more WSS variability in the descending aortic segments for the line-and voxel-to-voxel analysis, compared to the point analysis. Interestingly, these segments require less manual adaptation by the observer compared to the ascending aortic segments (at least for the distal descending aortic segment). The descending aorta (beyond the arch) is more fixed to the spinal column and therefore not sensitive to motion, as is the ascending aorta. For the descending aorta, we attributed the observed WSS variability to a higher velocity and more complex flow in these regions (compared to the ascending aortic segments), in the healthy volunteers.
These findings are important and give rise to reflection from a clinical perspective, for instance in patients with a stenotic bicuspid aortic valves or aortic dissections. Patients with stenotic bicuspid aortic valves often have ascending aorta dilatation/aneurysm formation [28]. Next to the difficulties for an accurate ascending aorta segmentation due to motion, higher velocities and complex flow phenomena are observed in the ascending aorta in these patients [29]. For the later mentioned, aortic dissection, the reconstructed geometry plays a key role in determining the complex blood flow in the true and false lumen [30]. Additionally, the presence of multiple intimal tears greatly influences the complexity of the flow in aortic dissections [31]. All of these characteristics lead to higher WSS variability depending on the applied WSS analysis method. Hence, it would be of clinical interest to perform an interobserver analysis in a subset of patients with a spectrum of aortic diseases (with and without high velocity outflow jets) with the proposed analysis methods, as a potential next step.

Limitations
A limitation of this study is a relatively small cohort of ten healthy volunteers with similar age. For a more robust analysis, a larger group of volunteers with wider age differences should be considered. Also, the study did not include patients with anomalies in the aorta, due to the ethical concerns with repeated examination. As has been shown, the lumen morphology has a great effect on WSS, which is often considered as one of the evaluation factors for certain arterial diseases [32]. Hence, for diseased patients, it might be crucial to segment the artery as close as possible to reality. Additionally, the diseased aorta may introduce more variability, not only by pathology, but also by the fact that 4D flow MRI acquisition is time consuming, which is always a more challenging examination in patients (more heart rate variability, less cooperation, difficulty not to move etc.) [10].
Several assumptions were made on the CFD part of this study, which may influence the results. First, is the assumption of laminar flow. Our operating range of Reynolds number was in the transition region and hence turbulence modeling should be considered and the effects of turbulence investigated.
Next, the rigid-wall assumption overestimates WSS, as demonstrated in Ref. [33]. However, this should not have an influence on our results, since the purpose of this study was to evaluate the WSS variability due to different segmentations in multiple volunteers. For our study, the actual movement of the aorta was similar, considering that the four compared aortas per each volunteer were based on two scans taken closely after each other and the physiological conditions (e.g. heartbeat) were very similar. Also, the segmentation was always performed at the peak systole. Additionally, the moving-wall approach is computationally very costly and requires additional patient-specific information, like the thickness of the wall and the elastic properties of the wall. All these parameters are difficult to obtain for each case and hence, using simulation approaches like Fluid-Structure Interaction or prescribed motion may lead to additional uncertainties in the simulations [33,34].
Third, we have considered just steady-state simulations with a blood flow simulated at the peak systole. The differences between steady simulations at the peak systole and genuine time-dependent simulations were addressed in Ref. [35]. Because of the steady-state simulations, some quantities used for evaluation of arterial flow cannot be obtained, e.g. time-averaged WSS and oscillatory shear index. However, since we have used the same boundary conditions for all simulations, the agreement in WSS between different segmentation should not be highly affected by the steady-state approach. Therefore, for this study the steady-state approach is sufficient.
Finally, a pre-processing procedure was applied for segmentations, which included cutting of the inlets and outlets, smoothing, and adding extensions. While we have kept all parameters the same and only one person was performing this procedure, some discrepancies can be still introduced to the surfaces due to the smoothing. For example, the positioning of the cutting planes on inlets could lead to a different size of inflow extensions. For a more rigorous analysis, the variability in preprocessing should be investigated as well.

Summary and conclusions
In the present study, we addressed the geometrically induced variability of the WSS in CFD-MRI coupled simulations. First, we adopted an approach often seen in literature, comparison of spatially averaged WSS mean in the five selected regions of the thoracic aorta (so-called 'Point Analysis' method). For both MRI and CFD results, agreement in global WSS for the different geometries is similar, showing, on average, good to excellent correlation in all selected parts of the thoracic aorta.
Next, we performed a more detailed visual and statistical analysis of CFD results. The circumferential averages of WSS were calculated alongside the centerline (so-called 'Line Analysis' method) and twodimensional mapping of the three-dimensional aorta wall values are performed to conduct voxel-to-voxel comparison (so-called 'Surface Analysis' method). In comparison to the 'Point Analysis', both 'Line' and 'Surface Analysis' show a lower agreement between different segmentations. The correlation varies between moderate and good for the 'Line Analysis', while it is between poor and good for the 'Surface Analysis' method. This reduced agreement is a consequence of minimal to no averaging in the 'Line Analysis' and 'Surface Analysis' approach respectively. Additionally, we observed lower correlation of WSS in the descending part of the thoracic aorta obtained from 'Surface Analysis' for various segmentations. This trend was not visible in the 'Point Analysis' and can be contributed to the more complex flow in this region.
Hence, our findings stress the importance of carefully analyzing the local WSS distributions of CFD simulations based on the 4D-Flow MRI segmentations. Finally, since we show that the WSS variability is similar for both rescan and intra/interobserver segmentations, CFD-MRI coupling shows the potential for studying the progression of aortic pathologies in serial follow up scans.

Funding
This research was financially supported by a grant from the Dutch Heart Foundation, The Netherlands (Grant Number CVON2018-08-RADAR).

Declaration of competing interest
None Declared. Contours of WSS for all four segmentations of one volunteer (V6) are shown in Fig. A.8 for both MRI and CFD. The scales for the contours were adjusted case-specifically: the WSS range is 0-30 Pa for the CFD results and 0-3 Pa for MRI results. Few similarities can be observed between MRI and CFD. WSS shows lower values in the ascending part of the aorta and gradually increases downstream for MRI as well as CFD. However, many differences can be found both in the localization of minima and maxima but also in the values. The maximal values of WSS are reaching just up to 3 Pa for the MRI, whereas for WSS based on the simulations, locations with WSS higher than 30 Pa can be found.

Appendix A.1Discussion
To establish the baseline of WSS CFD we compare the simulated results to WSSMRI for all four segmentations of one volunteer. Both methods lead to a similar global distribution of WSS with different absolute values (WSS CFD ∼ O (10) Pa and WSS MRI ∼ O (1) Pa). However, due to the assumptions in the simulation -rigid wall and no branching arteries -we are not able to perform a direct comparison. Both of these contribute to an overestimation of WSS CFD [36,37].
Moreover, it has been previously shown that the values of WSS MRI are underestimated [15]. To calculate the gradient of the velocity at the wall, the boundary layer has to be resolved adequately. The boundary layer thickness for the pulsating flows in an artery with diameter D is expressed as [22]: where δ is the boundary layer thickness. An order of magnitude estimate for the aorta gives the thickness of the boundary layer to be δ ∼ 0.1 mm. The resolution of the MRI measurements is Δx = 2.5 mm. Hence, MRI does not have a resolution high enough to adequately resolve the boundary layer and WSS cannot be properly estimated. Because of that, image-based WSS CFD should be always considered. However, for a patient-specific simulations, both the geometry as well as the boundary conditions should be considered as close to reality as possible. Study   Fig. B.9. Contours of WSS and the line for the data extraction for the three different meshes -Coarse (left), Normal (middle -used in all simulations) and Fine (right)for the geometry A of Volunteer 1 .   Fig. B.10. WSS alongside the out-seam in the descending thoracic aorta for three different meshes -coarse, normal and fine -for the geometry A of Volunteer 1

Appendix B. Grid Dependency
In order to perform a mesh dependency study three meshes were created for one of the cases: coarse, normal and fine. The coarse mesh, consists of 1.7 million elements, with the maximal cell size of 0.5D. The normal mesh consists of 4.2 million elements, with the maximal cell size of 0.3D. Finally, the fine mesh consists of 7.4 million elements, with the maximal cell size of 0.19D. For all meshes, the boundary layer settings were kept the same, as described in Section 2.2.2. The contours of WSS calculated on these three meshes are shown in Fig. B.9. No significant differences can be observed in contour distribution of the WSS at the thoracic artery wall. To asses in more details the local distributions of WSS, characteristic profiles were extracted along an arbitrary selected line along the descending part of the thoracic aorta (as indicated by black line), Fig. B.9. The extracted profiles of WSS are shown in Fig. B.10. Again, very small differences can be observed between the coarse and finer meshes. The maximal difference between the coarse and the normal mesh is 1.28% and between the coarse and the fine mesh 1.94%. The normal and the fine mesh agree well in WSS for most of the extracted data. The maximal difference between these two meshes is 0.64%. Based on all this, the normal mesh results were used for statistical analysis of the data.

Appendix C. Velocity information
Table C.7 shows the mean velocity and Reynolds number for the ten volunteers.