Improvement and validation of a computational model of flow in the swirling well cell culture model

Effects of fluid dynamics on cells are often studied by growing the cells on the base of cylindrical wells or dishes that are swirled on the horizontal platform of an orbital shaker. The swirling culture medium applies a shear stress to the cells that varies in magnitude and directionality from the center to the edge of the vessel. Computational fluid dynamics methods are used to simulate the flow and hence calculate shear stresses at the base of the well. The shear characteristics at each radial location are then compared with cell behavior at the same position. Previous simulations have generally ignored effects of surface tension and wetting, and results have only occasionally been experimentally validated. We investigated whether such idealized simulations are sufficiently accurate, examining a commonly‐used swirling well configuration. The breaking wave predicted by earlier simulations was not seen, and the edge‐to‐center difference in shear magnitude (but not directionality) almost disappeared, when surface tension and wetting were included. Optical measurements of fluid height and velocity agreed well only with the computational model that incorporated surface tension and wetting. These results demonstrate the importance of including accurate fluid properties in computational models of the swirling well method.

multidirectional flow that has recently been implicated in the localization of arterial disease (Mohamied et al., 2015). Furthermore, there are practical limitations to how long the shear can be applied; chronic application of flow, as occurs in vivo, is impracticable. To overcome this limitation, an increasing number of studies use the "swirling well" or "orbital shaker" method, in which cells are cultured on the base of a cylindrical dish or well while fluid motion is induced by placing the culture vessel on the horizontal platform of a shaker that translates in a circular orbit in the plane of the platform. The orbital motion forces a wave to swirl around the vessel; in broad terms, shear stress has lower magnitude and greater multidirectionality towards the center of the well than towards the edge. Local flow conditions are compared with local cell behavior. A recent review (Warboys et al., 2019) summarized 20 studies that have used this paradigm to investigate effects of shear on endothelial cells. Near the edge, cells demonstrate homeostatic properties such as elevated expression of eNOS and Krüppel-like factor 4 whereas cells towards the center exhibit potentially atherogenic behavior such as expression of intercellular adhesion molecule 1, vascular cell adhesion molecule 1, and E-selectin (reviewed in Warboys et al., 2019).
Because endothelial cells respond to the magnitude of the applied shear, the degree of its fluctuation, and changes in its direction, and because these factors have been implicated in triggering atherosclerosis, it is important to characterize the flow to which the cells are exposed with high accuracy. This is a more difficult problem for the swirling well than for the parallel-plate flow chamber and coneand-plate devices. Due to the complexity of the flow within the well and the presence of a free surface, the complete fluid behavior cannot be derived analytically, as with the other two devices, although recent advances permit the main flow features-the presence or absence of wave breaking and resonance-to be predicted for different combinations of orbital velocity, orbital radius, medium height, and well diameter. They also permit shear stresses at different locations on the base of the well to be approximated, but there is necessarily a compromise between tractability and accuracy (Alpresa et al., 2018a(Alpresa et al., , 2018b. Computational fluid dynamics (CFD) methods have therefore been used to model the flow and obtain the shear (Alpresa et al., 2018a;2018b;Bai et al., 2012;Barrett et al., 2010;Berson et al., 2008;Chakraborty et al., 2012;Discacciati et al., 2013;Filipovic et al., 2016;Ghim et al., 2017;Kim & Kizito, 2009;Salek et al., 2012;Thomas et al., 2017;Velasco et al., 2016;Warboys et al., 2010;Warboys et al., 2014;Zhang et al., 2009). Only a few of these studies have included surface tension and wetting in the model, and have experimentally validated the solutions.
Neglecting surface tension has been justified by computing the dimensionless Bond (Bo) and Weber (We) numbers, which indicate, respectively, the relative importance of gravitational acceleration to surface tension and the relative importance of inertia to surface tension: (1) We (2 ) , water 2 (2) where ρ is density, U is linear velocity, R is the radius of the well, σ is the surface tension coefficient, and g is gravitational acceleration (Salek et al., 2012). Values of both numbers have generally been >1.
Although ρ, σ, and g do not vary between different studies, U depends on the orbital shaker type and settings and R varies with the type of cell culture vessel. When the well diameter is large and the rotational rate and orbital radius are high, surface tension is negligible; gravitational and inertial forces largely determine the shape of the free surface. With a decrease in scale or linear velocity, surface tension plays a more important role. The contact angle between the liquid-air interface and the wall of the well also needs to be considered. It depends not only on surface tension but also on wettability. For a water-air interface in contact with a hydrophobic material, angles are >90°, and for a hydrophilic material they are <90°, giving rise to concave and convex menisci, respectively.
Here we used CFD with and without surface tension plus wetting to model fluid flow in a configuration that has been used in several published studies (Arshad et al., 2021;Ghim et al., 2017Ghim et al., , 2021Potter et al., 2011Potter et al., , 2012Warboys et al., 2010): cells cultured in a standard 12-well plate (or on a Transwell ® filter insert in a six-well plate) on a shaker with an orbital radius of 5 mm and rotational rate of 150 rpm.
Incorporating surface tension and wetting changed the computed flow regime from a breaking to a non-breaking wave and virtually abolished the center-to-edge change in shear magnitude, but not shear directionality. The results were validated by experimental Both were achieved optically, the former by measuring the absorbance of light passing vertically though the fluid and the latter by PIV.
Numerical characterization involved models of the swirling well that differed only in whether they did or did not include surface tension plus wetting. These physical properties of the fluid were determined experimentally, using a tensiometer to assess surface tension and Meniscus heights were corrected for refractive index errors using a calibration factor obtained by comparing known heights with heights measured by the confocal. To compute the contact angle, a theoretical meniscus height y at radial location x from one wall of the well, relative to the height at which the liquid pressure is equal to atmospheric pressure, was fitted to the experimentally determined shape. The theoretical curve was the sum of two exponentials: and ρ is the density, S is the surface tension, g is the gravitational acceleration, θ is the contact angle, and R is the radius of the well. θ was varied until a best fit was achieved. A plot of this equation is given in Figure S1. The two exponentials in Equation 3 essentially give the two sides of the meniscus, one decreasing in height with increasing x and the other increasing, and contributing equally in the center. The precise curvature is defined by ρ, S, g, and θ, which are incorporated into the two constants, C 1 and C 2 .

| CFD modeling
The motion of fluid in one well of a 12-well plate translating in a circular orbit on a horizontal platform was simulated using commercial CFD software. The geometry consisted of a cylinder having diameter 22.1 mm, height 10 mm and an open surface, filled with liquid to a mean depth, d, of 2 mm ( Figure S2). The remaining space was filled with air. The base and side walls of the cylinder were defined as no-slip boundary conditions and the open surface was defined as a pressure outlet (atmospheric pressure). The motion of fluid in the well was described by continuity and Navier-Stokes equations: where u is the velocity, ρ is the fluid density, t is the time, p is the pressure, ν is the kinematic viscosity, and f is the force acting on the fluid. The equations assume that the fluid is incompressible and Newtonian. They were converted to discrete finite difference equations and solved using STAR-CCM+ (version 11.02.010-r8).
An acceleration given by: was applied in the x, y, and z planes, respectively, with −9.81 being the gravitational acceleration, A the orbital radius (5 mm) and ω the angular velocity (15.7 rad/s).
A volume fraction, C, of the reference phase (liquid) was specified for each grid cell. Grid cells filled with liquid had C = 1; those filled with air will had C = 0, and cells containing the interface had values in the range 0 < C < 1. The free surface was defined as the isosurface where C = 0.5. Velocity and pressure were obtained using the segregated flow model, which solves the momentum and continuity equations sequentially in an uncoupled manner, linking them using a predictor-corrector approach. Physical properties of the various liquids incorporated into the computational models are given in  (Alpresa et al., 2018a).
The "multiphase interaction model" was used to incorporate surface tension. Following Brackbill et al. (1992), the surface tension force (f σ ) acted normal to the free surface (n) and depended on the surface tension coefficient (σ) and the mean curvature of the free surface (K): This led to artefacts: parasitic currents at the interface gave rise to velocities that were much larger than those in the rest of the fluid.
They were damped using an Interface Momentum Dissipation (IMD) model, which adds artificial viscosity at the surface. Large values of IMD can have a detrimental effect on wave behavior. There is no general cut-off value for IMD; it needs to be determined on a caseby-case basis. Figure S3 shows that the addition of 1.0 IMD damping to the CFD model virtually eliminated the excessive mean and average velocities seen at the surface with lower or no damping.
A value of 3.0 was chosen; it had no detrimental influence on the overall pattern of flow but improved convergence of the solution compared to a value of 1.0.
The IMD value of 3.0 was used in conjunction with a surface tension of 72 mN/m and a contact angle of 70°, which is the contact angle of water with poly(methyl methacrylate) (PMMA). The time step was reduced to 5 × 10 −5 to satisfy the CFL criterion. The number of mesh cells was maintained at 652,000 following a further mesh convergence study which showed that it gave mean WSS values indistinguishable from those obtained with 1,350,000 cells and maximum WSS values that were only around 10% lower than with the finer grid ( Figure S4).
In additional simulations, the model with surface tension and wetting was run with the rpm of the shaker altered from 150 to 100 rpm, and with the fluid properties changed from those of water to those of ethanol (Table 1). The contact angle for ethanol is 20°. were extracted, the liquid domain was sliced into 0.05-mm-thick layers using ParaView (v.5.4.1; www.paraview.org), and the velocity information averaged.

| Synthetic particle image generation
Using a custom Python script, the liquid domain of the CFD model was randomly seeded with 200,000 or 500,000 particles. The particle positions were then displaced linearly according to the velocity field using ParaView. The positions of the original and displaced points were processed through the PIV Synthetic Image Generator (SIG) (Lecordier & Westerweel, 2004) to produce an image pair. The settings for the SIG were kept the same as the experimental set-up.
A range of different parameters were tested using the SIG. The thickness of the laser sheet was varied between 3.75, 7, and 10 mm, the laser sheet profile was Gaussian or uniform, and the location of the laser sheet was moved +1, +2, or −1 mm from its original position.
All of these parameters were tested while varying the time between laser pulses from 0 to 4000 µs. Figure S5 illustrates the PIV apparatus. An optical dovetail rail (Thorlabs) was screwed onto the platform of the orbital shaker. It was used to secure a high-speed camera (Dantec Dynamics FlowSense EO) and a custom-made PMMA well sitting on a 3D printed frame that held a mirror at 45°to its base. The PMMA well had the same dimensions as one well of a Corning Costar ® 12-well plate. A 532-nm laser (Laser Nano L laser head with an LPU550 PSU power head), synchronized with the camera and having a pulse interval of 1 ms, was placed 1 m away. A silica plano-concave cylindrical lens with a focal distance of 500 mm and a silica plano-convex lens with a focal distance of −25 mm were used to a produce a flat laser sheet parallel to the base of the well.

| Experimental PIV setup
To overcome reflections, the back inner wall of the well was sprayed with red-orange matt fluorescent paint (Rust-Oleum) which fluoresced at ≥600 nm. A 532 ± 2 nm band-pass filter was placed in front of the camera | 75 lens to pass the 532-nm laser light scattered by the particles while rejecting the higher wavelength fluorescence from the paint.
The water in the well was seeded with polyamide particles having a mean diameter of 5 µm, diameter range 1-10 µm, and density 1003 kg/ m 3 (Dantec Dynamics). The hydrophilicity of the particles was increased by plasma treatment for 1 min before dispersion.
To ensure that the wave was captured in the same position in each orbit, and hence to allow averaging of images, an external trigger signal was generated using a metal strip secured to the shaker platform so that it crossed an Arduino photo interrupter module once every rotation. The maximum camera trigger rate was less than the shaker rotation frequency, so every third rotation was captured.

| PIV processing
Image pairs were post-processed using PIVlab, a Particle Image Velocimetry tool for MATLAB (Thielicke & Stamhuis, 2014). A circular mask was used to isolate the area of interest and contrast-limited adaptive histogram equalization (CLAHE) was applied to improve the image quality.
Each image was sub-divided into interrogation windows of 64 × 64 pixels (pixel size = 6 µm). Two passes were used, 64 × 64 and 32 × 32 with a step size of 50%. Interrogation windows were cross-correlated between image pairs by fast Fourier transform window deformation: where i and j are the length and width of interrogation window and C is the correlation matrix giving the best displacement of particles between interrogation windows A and B. The displacement between an interrogation window in the first frame and the window from the second frame giving the highest cross-correlation with it, divided by the time delay between the images (which is equal to the time delay between laser pulses), gave a velocity for each window. Velocity outliers (magnitudes > 0.22 m/s) were replaced by interpolated vectors. One image pair was used to produce each 2D vector map for the experimental data and an average of 10 pairs for the synthetic data.
As with the CFD-ST and CFD+ST models, PIV experiments were also conducted after the rpm of the shaker was changed from 150 to 100 rpm and after the liquid was changed from water to ethanol.

| Measurement of fluid height in a swirling well
To map the height of the fluid in a swirling well, 14 mg/ml EBD was added to the liquid in the well, a light sheet was placed under the well to trans-illuminate the fluid from below and a camera (Samsung S8) was secured on the orbital shaker platform 10 mm above the center of the well ( Figure S6). The camera was focused on the bottom of the well. Imaging parameters were: sensitivity, 100 ISO; aperture, f/1.7; white balance, 5400 K; exposure time, 1/180 s; image size, 4032 × 1969 pixels.
Initially, the well was stationary and images were taken with decreasing volumes of EBD solution, from 2304 to 192 µl in increments of 192 µl (corresponding to a height of 6-0.5 mm in increments of 0.5 mm, in the absence of a meniscus). Images were then taken of the swirling well with the platform translating at 100 rpm or 150 rpm, using 767 µl of EBD solution (corresponding to a mean height of 2 mm).
Images were post-processed in MATLAB R2016a. Pixel intensities from the images of the static well were inverted and fitted

| Measurements of surface tension and contact angle
Surface tension data are shown in Figure S7.    comparison.
An increase from 200,000 to 500,000 particles produced negligible differences in the average velocity but dramatically reduced the large errors in maximum velocity seen at the lower seeding density for pulse intervals below 1000 µs.
Changing the light sheet position and profile had relatively small effects; in general, the highest errors were seen at z = +1 mm. (This F I G U R E 4 Effect on average and maximum velocity, obtained in the synthetic PIV solutions, of changing the properties of the laser sheet (vertical location, thickness, and beam profile), the interval between laser pulses, and the particle seeding density. Unless otherwise stated, the laser profile was Gaussian and the sheet was 3.75 mm thick. Horizontal black lines indicate the mean and maximum velocities (0.19 and 0.063 m/ s, respectively) obtained directly from the corresponding CFD solution, which did not included surface tension and wetting. Vertical dashed black lines indicate the laser pulse interval chosen for experimental PIV. CFD, computational fluid dynamics; PIV, particle image velocimetry position refers to the bottom of the light sheet). Switching to a thicker sheet (7 instead of 3.75 mm) had a negligible effect. Using a uniform rather than a Gaussian beam profile did have an effect, but only for the maximum velocity when detected with the high seeding density and a delay of around 1000 µs, where it increased the discrepancy with the CFD result. Thus a laser sheet with Gaussian profile and thickness of 3.75 mm appears sufficient despite the maximum height of the liquid reaching 4.2 mm in the computational models, and the results are insensitive to the precise position of the sheet so long as the first millimeter above the base of the well is included.
The interval between pulses had by far the largest effect. The average velocity was highest at 10 µs but dropped by 100 µs, beyond which it remained fairly constant and close to the CFD value, until 1000 µs; it continued to drop as the time between laser pulses was increased further. The maximum velocity was also closest to the maximum CFD velocity at ≈1000 µs but increased sharply prior to 1000 µs as well as for longer delays.
Thus 1000 µs seem to be the optimum interval for accurately capturing the flow fields. At this delay, there was negligible effect of seeding density on either the average or the maximum velocity and, as already noted, the height of the laser sheet also had only minor effects with the exception of position z = +1 mm.
Regardless of the number of particles and the position or profile of the laser sheet, the maximum velocity at 1000 µs was always slightly higher than the result obtained directly from the CFD model.   Figure 6. Ethanol was used to reduce surface tension (and is also likely to have affected particle dispersion) and 100 rpm was used to eliminate the breaking wave even in the CFD-ST model. Under both conditions, there appeared to be a better agreement between the experimental PIV and the CFD+ST model than was observed for the baseline (150 rpm, water) case shown in Figure 5.
(Note that there is a readily apparent discrepancy between CFD-ST and CFD+ST simulations for ethanol in Figure 6, as there was for the baseline case shown in Figure 5).
The average and maximum velocity magnitudes calculated for each map are given in Table 2 to that simulated with surface tension and wetting ( Figure 6). In neither case is the disagreement as large as for the baseline model, presumably reflecting the fact that the arc of rapid flow at the top of the well is seen in both CFD+ST and PIV in the ethanol case, and in neither for the 100 rpm case.
The experimental PIV maps were rotated by 0-360°and the correlation coefficient with the corresponding CFD models was calculated. Separate Pearson's coefficients were calculated for the velocity magnitude and the vector orientation since the two gave best matches at different rotations; the average of the two correlation coefficients is shown in Table 3. For the baseline (150 rpm, water) and ethanol simulations, experimental PIV correlated better with the CFD +ST than the CFD-ST model. (The coefficients were approximately 50% higher). Both correlation coefficients for the ethanol simulations were higher than both coefficients for the baseline simulations.
Experimental PIV correlated better with the CFD-ST than the CFD+ST model in the 100 rpm simulations; nevertheless, both coefficients were still higher than both coefficients for the baseline case.

| Experimental versus computational height maps
The calibration curve for the relation between liquid height versus pixel intensity, obtained by microscopy of static wells containing various volumes of a dye solution, is shown in Figure S8, together with the exponential curve fitted to the data. The exponential was truncated at 10.4 mm from the edge of the well to avoid errors caused by residual reflections from the side wall.
The calibration was then applied to obtain the height of the liquid across wells that were swirled at 150 or 100 rpm, and the results were compared to corresponding CFD-ST and CFD+ST height maps ( Figure 7). To match the phases between computed and experimental height maps, the experimental height maps were rotated as above Note: The maps were rotated to achieve the best correlation coefficient for either the vector magnitude or the vector direction, and the average of these two r values is tabulated.
F I G U R E 7 (a) Height of the free surface obtained by experiment and by CFD simulations with and without surface tension plus wetting (CFD+ST, CFD-ST), and corresponding pixel-by-pixel scatter plots and their correlation coefficients for (a) baseline and (b) 100 rpm cases. The black circles are drawn at a radial distance of 10.4 mm from the center of the well; experimental data beyond that are less accurate due to optical reflections from the side wall. Note that wave breaking, evident as a disturbed leading edge to the wave, is only apparent in the baseline model. CFD, computational fluid dynamics; PIV, particle image velocimetry and the rotation that gave the highest Pearson's correlation coefficient between the two was used.
At both rotational speeds, the experimental data agreed better with the CFD+ST height map than the CFD-ST one. (The correlation coefficients are also shown in Figure 7). The correlation coefficients were lower at 100 than 150 rpm. The overall correlation coefficients at 150 rpm were >0.9 with or without surface tension plus wetting.
Consistent with the velocity maps, the maps from CFD-ST simulations showed less curvature than the CFD+ST maps or the experi-  Figure 7a). CFD, computational fluid dynamics; PIV, particle image velocimetry numerically simulated a wide range of swirling well configurations and found that the wave would break if: where  fig. 11b and d), as predicted by the criterion given in Equation (8) and by the modified equation given at the end of the Conclusions section, below. A more recent study carried out by Thomas et al. (2017) compared CFD data with PIV measurements that spanned the entire diameter of the well. Velocity vector components over one cycle differed <5% between CFD and PIV, and magnitudes averaged over the inner 20% of the radius differed by only 2.4%. However, the mean fluid depth in the experimental study was 10 mm, which is fivefold greater than in our study; Equation (8) suggests that increased mean fluid height reduces the likelihood of wave breaking.
In our study, wave breaking was clearly visible in the height maps produced by a conventional CFD model and was eliminated when surface tension and wetting were incorporated. The modified model also gave greater height near the edge of the well and moved the minimum height away from the wall, and mean and maximum velocities within the fluid were reduced except in the boundary layer at the base of the well.
Modifying the model had a major impact on WSS at the base of the well. Its magnitude was reduced and its pattern at the wavefront was altered. There was negligible difference in the plot of mean WSS vector direction versus radial distance, but differences were seen in the instantaneous vectors: flow was purely multidirectional close to the center but had greater magnitude when surface tension and wetting were included, while towards the edge the modified model gave more uniaxial behavior and greater reverse flow. There were corresponding changes in the dependence of various WSS metrics on radial location. The relative changes in TAWSS and transWSS are of particular significance. If the well is divided into edge and center regions having equal area, which can be achieved by using a radial distance of 7.8 mm as the boundary, then the difference in TAWSS between the two regions is reduced from being 30% lower in the center to being only 7% lower, whereas transWSS is hardly changed (58% higher at the center without surface tension and wetting, and 53% with them). Differences in cell behavior between the two regions observed in previous studies are thus more likely to reflect changes in multidirectionality than in shear magnitude.
Before testing the CFD results against experimental PIV, we first optimized the parameters to be used in the PIV experiments by using "synthetic PIV," a computational model of PIV generated by applying the synthetic image generating software of Lecordier and Westerweel (2004) to our CFD simulations. Comparing the mean and maximum velocities obtained by synthetic PIV with those obtained directly from the underlying CFD showed that the critical variable was the interval between laser pulses. When that was 1 ms, agreement was good and relatively insensitive to particle density or the vertical position, thickness, and profile of the laser sheet. (Failing to capture flow near the base of the well was the next most critical issue).
Using optimum values of these parameters, the synthetic PIV and the CFD simulations from which they were generated gave essentially identical results, both when surface tension and wetting were incorporated and when they were not. Even with velocity changing almost fivefold across the well, the postprocessing algorithm was able to compute accurate flow fields using a constant interrogation window size and time between laser pulses. The maps of velocity vectors, averaged through the thickness of the laser sheet, were also consistent with the WSS maps described above and, like them, showed substantial differences between the models with and without surface tension plus wetting: the magnitudes were lower in the modified model, there was a different pattern of vector direction and magnitude at the front of the wave, and the positions of the source and sink of the streamlines were changed.
The good agreement between data obtained by synthetic PIV and directly from the underlying CFD leads to the expectation that experimental PIV should provide a true record of the real flow in the swirling well, subject only to experimental errors. A key finding was that the agreement between experimental PIV and CFD was substantially better when surface tension and wetting were added to the baseline model. Nevertheless, discrepancies remained: the experiments gave lower overall velocity magnitudes, there was a shift of an arc of high magnitudes away from the edge of the well, and the source of streamlines was further from the edge of the well.
When ethanol was used instead of water, the correlation coefficients were nearly doubled. The improved correlation may indicate an error in how the CFD handles surface tension, since ethanol has a lower value than water, which was used in the baseline configuration.
However, the correlation remained higher for the model with surface tension plus wetting than for the model without them. An alternative explanation is that the seeded particles did not mix well with water despite plasma treatment, leading to clumping; this was occasionally visible by eye. Particles would disperse better in ethanol. Larger particles would hinder flow tracking (Melling, 1997) and could also lead to alteration of the flow field.
When the rotational speed of the orbital shaker was reduced from 150 to 100 rpm, the map of velocity vectors was markedly of 100 rpm and an orbital radius of 6 mm; these values are reasonably close to those studied here. The solutions showed "splashing" for certain fluid heights, but the simulations appear only to consider the first quarter of an orbit after starting the shaker, when it is unlikely that a steady-state has been reached. A numerical model was developed by Toro et al. (2018), based on equations from Lubarda (2013), to derive free surface profiles, and was applied to 96-well and 384-well microplates on an orbital shaker rotating at rates >500 rpm.
When surface tension was included, there was good agreement with surface profiles obtained by imaging over a wide range of parameter values. There was poor agreement in the absence of surface tension.
Furthermore, numerical solutions for dimethyl sulfoxide, which has a surface tension of 48 mN/m, gave surface profiles that were markedly different to those for water (72 mN/m), again suggesting that surface tension has a large effect. Toro et al. (2018) also modeled two contact angles-87.4°and 102.1°-to represent water in, respectively, polystyrene and polypropylene wells. This had a much smaller influence on surface profiles: for example, surface tension had a >20-fold larger influence on maximum wave height (>200% change vs. <10% change). Surface tension acts to minimize the free surface area of a liquid and, for swirling liquid in a cylindrical container on an orbital shaker, will therefore act to flatten the surface. The contact angle governs the profile at the point where the free surface intersects the side wall.
Changes in contact angle can influence the profile further in, but are unlikely to have the same dampening effect on wave behavior.
The present study suffered from a number of limitations: 1. Experimental height measurements agreed better than experimental PIV data with the results of CFD+ST simulations. This is consistent with the presence of experimental errors in the PIV.
We have discussed above the possibility that the seeded particles clumped together. Other potential sources of error are the thickness of the laser sheet and the failure to capture velocities in the z-plane, although these limitations were also present in the synthetic PIV, which agreed well with the numerical results. 4. Surface tension and contact angle were not changed individually and hence their effects cannot be separated. However, as noted above, Toro et al. (2018) showed that influences of surface tension on the surface profile were >20-fold larger, and we would expect that discrepancy to be even greater in the present study because it used bigger wells in which less of the fluid behavior would be affected by the walls. In particular, the reduction in wave breaking that we observed over large portions of the well diameter in the modified model is plausibly explained by the dampening effect of surface tension.

| CONCLUSIONS
Surface tension and wetting should be included in simulations, despite the >10-fold increase in computational cost (>10 days vs. <1 day on the Imperial College HPC system), unless they are already known to have only negligible effects on the parameter values of interest. The Weber and Bond numbers do not appear to be a reliable indicator of when this is likely to be the case.
The baseline parameter values used in the present study have previously been employed in earlier investigations of cell behavior, for which WSS was computed using models that did not incorporate surface tension and wetting. The different results obtained here with the modified model alter the conclusions of those studies. In particular, they increase the likelihood that radial variation in cell behavior reflects changes in transWSS rather than TAWSS. suggest that the threshold value needs to be increased; breaking was not seen even at Δh/(2 R) ≈ Γ.