Immersed boundary method applied to flow past a tree skeleton

.


Introduction
Small-scale tree configurations, such as solitary trees, alleys, wind breaks, and small forests occur frequently in both urban and rural landscapes.Their effect on the immediate environment is significant, both in terms of wind, temperature and moisture level.Therefore, the impact of trees is important in the context of, for example, urban comfort (Hefny Salim et al., 2015;Matsuda et al., 2018) and agriculture (Wang and Takle, 1996).During severe wind conditions, trees can also pose a severe danger, since breakage can cause material damage and death.For these reasons, it is important to accurately represent the effect of trees in wind and meteorological models.This is challenging because of their complex multi-scale structure (Gosselin, 2019), which requires proper treatment of a wide range of length scales.
An early attempt to numerically simulate the effect of an idealized solitary tree was made by assuming an aerodynamic porosity and a tree geometry approximated by a sphere or a cone (Gross, 1987).However, the most common way of representing a single tree in wind models is via a drag force parameterization (Dellwik et al., 2019;Endalew et al., 2009;Poh et al., 2020), with inspiration from homogeneous (Finnigan, 2000) and heterogeneous (Boudreault et al., 2015;Schlegel et al., 2012) forest canopy parameterizations.With this parameterization, the tree geometry is spatially averaged and represented with an area a per unit volume.This area is then multiplied with a drag coefficient in the common drag equation.The drag force parameterization mimics the flow through the porous tree crowns and foliage in that it controls the resistance of the flow through it.Although this approach can represent the flow past trees quite accurately, the dependence on the drag coefficient is a drawback.In contrast to the area a, the drag coefficient is not a directly measurable characteristic of the tree and, therefore, has to be empirically determined either by calibrating against a measured drag force on the tree (Dellwik et al., 2019) or by calibrating the drag coefficient to match observed wind fields (Endalew et al., 2009;Wilson and Shaw, 1977).Since drag force observations on trees are rare, and furthermore may show significant variability between trees of similar stature (Manickathan et al., 2018;Rudnicki et al., 2004), this approach may lack quantitative accuracy.
A more direct way of estimating the drag force on the tree is to resolve the structure and insert hard impenetrable surfaces into the numerical solver.Due to computational limitations, this approach is currently only practicable for the main part of the tree structure, such as the stem and the major branches.This approach was used by (Endalew et al., 2009) who made RANS simulations on two scaled artificial plastic trees where the stem and branches were resolved with an unstructured grid while the foliage was modelled by introducing drag forces in porous sub domains surrounding the branches.They compared the results of this integrated method (IM) to those obtained when modelling the entire tree with forcing terms and found that the IM more accurately reproduced the wake characteristics measured behind the trees in a wind tunnel.A challenge with this approach is the grid generation, where surface irregularities and intersection of scales in the complex tree structure may challenge the capacity of automated solutions.
A more flexible and direct representation of the irregular stem and branches can be provided by the immersed boundary method (IBM) (Peskin, 1972) which does not require the grid to conform with the simulated surface.This approach was adopted by (Chester et al., 2007) and further developed by (Graham and Meneveau, 2012) in their large simulations of the flow past fractal trees.They proposed a so-called renormalized numerical simulation (RNS), which elegantly splits the forces into a part due to the large scales of the tree and a part due to the small branches.The IBM then directly resolves the large scales, while the latter unresolved part is modelled using only information about the forces on the resolved part.In this way, the method is consistent with the notion of LES where information from the large resolved scales of the flow are used to determine the subgrid stresses.Besides showing that the drag predicted by the RNS is not very sensitive to a range of tested grid resolutions they also found it to be in rather good agreement with the drag measured in a water tunnel.Whereas this result points to how the smaller scales can correctly be parameterized for fractal structures, it is less clear how suitable the method is for specific real trees, that may not show typical fractal behavior.
As a compromise between the surface-conforming grid approach by (Endalew et al., 2009) on a real tree, and the IBM method by (Chester et al., 2007) and (Graham and Meneveau, 2012) on tree-like fractal structures, we here simulate a tree skeleton deduced from a high-detail scanning of a real tree.We use IBM together with both Reynolds averaged Navier-Stokes (RANS) and detached eddy simulations (DES) of the flow field.We focus the RANS simulations only on predicting the forces on the tree, while the DES also is aimed at simulating the turbulent wake.The tree skeleton was placed in a large scale wind facility, and measurements of both the drag force and the wind field in its wake were recorded.Similar to (Chester et al., 2007), the near-surface flow is parameterized via a surface roughness, and we can adjust this value to fit with the observed drag force at a given wind speed, such that we ensure that the tree exerts the correct overall drag force on the flow.The main aim of the study is to evaluate the performance of the IBM for a complex, non-fractal tree skeleton structure in a turbulent inflow.However, since the surface roughness turned out to be more important than expected, the sensitivity to this parameter is also investigated by means of numerical simulations.

Experimental method
The experimental part of the present work was conducted in the WindEEE Dome (Hangan, 2014) where measurements of both the forces on and flow field around a small scale tree were conducted at different inflow conditions.In the following the different parts of the experimental campaigns will be described.The description of the used instruments will be brief but more details about them can be found in the work by (Enus ¸ et al., 2020).

Tree model
The tree used in the experiment has a height of H = 0.41 m and is a simplified and down-scaled representation of the 6.4 m tall tree which was the subject for the full scale wake investigations by (Dellwik et al., 2019).The impact of the down-scaling is studied numerically in appendix and it is shown that there are no significant scaling issues associated with the results obtained on the model tree.The process of generating the scaled tree involves scanning the real tree with a high resolution lidar as described by (Dellwik et al., 2019) and then converting these scans into a watertight model surface representing the stem and main branches using the method described by (Gabriel, 2017).Finally, the model tree was 3D printed in PLA plastic and placed in the center of the WindEEE Dome as shown in Fig. 1.Note that due to difficulties in the 3D printing process, the model tree is characterized by significant surface imperfections which can clearly be seen in the right plot of Fig. 1.As discussed in Section 5, these imperfections have a big impact when parameterizing the surface roughness used in the numerical model.

Flow measurements
Measurements of both the inflow and the wake were conducted at different wind conditions by regulating the wind tunnel fan speed and a realistic turbulent inflow was generated by introducing 0.15 m high roughness elements in the region upstream of the tree (see Fig. 1, left).
The measurements were recorded using a rake of four-hole Pitot tubes (Cobra Probe 100 series, Turbulent Flow Instrumentation Pty Ltd., Australia), which can measure wind speeds between 2 m /s and 100 m /s with an accuracy of 0.3 m/s and flow at angles up to 45 ∘ with an accuracy of 1 ∘ .When the fan speed was running at 100%, the measured wind speed generally never went outside the range allowed by the Pitot tubes.However, when the fan speed was lowered (especially below 60% of it's maximum), then the velocity in the wake more often went below the lower limit of the Pitot tubes and therefore it was decided only to use the measurements obtained at a fan speed of 100% in the analysis of the wake.The measurement campaigns at lower fan speeds could however still be used in the analysis of the forces as described in the Section 2.3.The wake was measured at cross-sections located x/H = 0.63 and x/H = 1.90 downstream of the tree as sketched in Fig. 2. Each cross-section spanned an area defined by − 0.49 ≤ y/H ≤ 0.49 and 0 ≤ z/H ≤ 1.95, where y and z denote the lateral and vertical direction, respectively and the tree is located at (x, y, z) = (0, 0, 0).The inflow conditions were estimated in two different ways: In the first approach, measurements were conducted in an empty tunnel at the exact position where the tree would later be mounted.The second approach used the wake measurements obtained at y/H ± 0.49 since at these positions, the tree was expected to have only a minor effect on the flow (something which is confirmed later in Figs. 10 and 11).An analysis of the inflow measurements revealed that the flow in the tunnel is not typical for a standard neutral atmospheric boundary layer, but rather reflected a rough to smooth land surface transition, with less shear and turbulence in the lower part of the measurement range, and higher above.This is Fig. 3. Measured mean wind speed (left) and standard deviation (right) of the inflow normalized with u 0 , where u 0 denotes the free-stream velocity at z /H = 0.74.The shaded blue area indicates the vertical height range of the lowest to the highest branch.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)illustrated in Fig. 3 which shows the mean and standard deviation of the undisturbed velocity measured at the different positions when the fan is operating at 100%.The mean velocity profile is far from logarithmic and the turbulence intensity varies significantly with height.The imbalance in the inflow is clearly seen and leads to an acceleration of the flow of approximately 10% from x/H = 0 to x/H = 1.90 and a turbulence decay of a similar order of magnitude (Fig. 3, right).This lack of inflow equilibrium is a major disadvantage of the experiment and reduces the level of agreement that we can expect when validating the simulations.However, since the flow undergoes only minor development from x /H = 0.63 to 1.90, it was decided to use these downstream measurements as representative for the "true" inflow when comparing with simulations.

Force measurements
The forces on the tree were measured concurrently with the wake measurements using a force moment sensor (JR3 45E, Woodland, Ontario Canada) on which the tree was mounted.The balance measures the forces and moments in three directions but here we will only consider the force in the streamwise direction because it is the dominant component.The sensor samples at 1000Hz and can measure streamwise forces in the range ±245N with an accuracy of 0.12N.

Numerical method
The flow field is simulated using both steady state Reynolds Averaged Navier-Stokes (RANS) and Detached Eddy Simulations (DES) while the tree is represented using an immersed boundary method (IBM).The idea of using RANS is to solve directly for the steady state solution, whereas DES is an unsteady model, which combines RANS near the tree surface with large eddy simulation (LES) of the flow away from the tree.The advantage of using RANS is that it is much faster than DES.However, whereas we find that the two models predict similar forces on the tree, it is well known that RANS is often inadequate for representing wakes in the atmospheric boundary layer as shown by (Laan et al., 2015).Therefore, we will in the following only use RANS for predicting the drag of the tree, whereas we will use DES for the numerical predictions of the wake.It should mentioned that with the setup used in the following, the DES is effectively LES everywhere except for the first few cells off the tree surface which is modelled as RANS.In the following we will describe the different parts of the model setup.

Flow solver
The simulations are carried out using the Navier-Stokes solver EllipSys3D (Michelsen, 1992;1994;Sørensen, 1995), which approximates the incompressible RANS using a finite volume discretization on block structured collocated grids.We solve the coupled momentum and pressure correction equations using the semi implicit method for pressure linked equations (SIMPLE) algorithm (Patankar and Spalding, 1972) combined with the modified Rhie/Chow interpolation method proposed by (Troldborg et al., 2015) to avoid pressure-velocity decoupling.
As mentioned above we use two different turbulence models: For the RANS we use the k − ω shear stress transport (SST) model by (Menter, 1994) and for the DES we use the version proposed by (Strelets, 2001).
For the RANS simulations we discretize the convective terms using the quadratic upstream interpolation for convective kinematics (QUICK) scheme (Leonard, 1979).For the DES, we use a hybrid scheme, which switches between the QUICK scheme in the RANS regions and a fourth order CDS scheme in the large eddy simulation (LES) regions.The unsteady DES solution is advanced in time using a second order three-level implicit time stepping method which relies on sub-iterations at each individual time step.The time step used for the DES is Δt = 2.5⋅10 − 4 s.

Immersed boundary method
In the present work we use a feedback forcing approach (Goldstein et al., 1993;Saiki and Biringen, 1996) to determine the force exerted by the immersed boundary (IB) on the flow field.Thus, the boundary force is defined as a singular function along the boundary which is then projected to the background computational grid via a certain regularized delta function.The transfer function can be defined in many ways but here we use the modified Rhie/Chow interpolation method by (Troldborg et al., 2015), which ensures a consistent treatment of the boundary force in both the momentum and pressure correction equations.
Our implementation of the immersed boundary method follows approximately the same steps as shown in previous work (Iaccarino and Verzicco, 2003), however, it has some unique features and therefore we will describe the implementation in more detail below.Besides we will also comment on some of the known drawbacks of the method as discussed by (Iaccarino and Verzicco, 2003) and how we in part can overcome them.Although we find that the method is not problem-free, it turns out that it works rather well for simulating the tree considered in the present work.

Representation of the immersed boundary
The surface of the immersed boundary is discretized into a number of triangular shape cells (SC) which together form a watertight surface.Each triangle has an outward pointing unit normal vector n SC as sketched in Fig. 4(a).
The domain cells (DC) in the computational grid that are intersected by the IB surface are flagged as forcing cells (FC), i.e. these are the cells in which a force term is added in order to fulfill certain boundary conditions (See Fig. 4) As mentioned above, the force from the IB is projected onto the computational grid by using a modified Rhie/Chow interpolation technique (Troldborg et al., 2015).In this connection we need to determine the point of action of the force inside each FC.We call this point the forcing point (FP) and define it as the center of gravity of the part of the IB which intersects the FC.Besides the FP we also need to compute S FC , which is the area of the part of the IB surface located inside the FC and n FC which denotes its outward pointing unit normal vector.The area S FC is needed because the force applied in each FC depends on how much of the IB is inside of the FC, while n FC is used in connection with turbulence modelling and near wall treatment.Finally, a user defined probe point (PP) is also defined for each FP, which is located at a distance of d PP from the FP on a line perpendicular to the IB (see Fig. 4c).The PP is used for setting boundary conditions for the shear stress and turbulence quantities as described in Section 3.2.2.

Calculation of the boundary forces
The force exerted by the IB on the flow is decomposed into components which are normal and parallel to the IB.The normal force component is determined using a feedback forcing approach (Goldstein et al., 1993;Saiki and Biringen, 1996) where the idea is to compute forcing terms in response to the difference between the actual and desired velocity at the IB and insert these in the governing flow equations and thereby drive the velocity through the IB asymptotically to zero.In practice the forcing is computed via a PI control strategy: Here K P [kgs − 1 ] and K I [kgs − 1 ] are constants, subscript i denote the iteration number, Δ is the number of iterations used for evaluating the integral term and e i = v FP n,i − v n,0 is the difference between the normal velocity at the forcing point and its desired value v n,0 .Note that v FP n is determined using linear interpolation from the surrounding cells and that v n,0 = 0 because the present work considers only non-moving boundaries.In the following we use Δ = 20 for both steady and unsteady simulations.Since the unsteady simulations are conducted with 5 sub-iterations per time step, this is equivalent to evaluating the integral term over 4 time steps.
The force f FP n acts at the forcing point and we need a method to transfer this force to the computational grid.As mentioned above we use the modified Rhie/Chow method by (Troldborg et al., 2015), which ensures a strong coupling of the velocity and pressure in the vicinity of the IB.
The dependence on the coefficients K P and K I is one of the biggest disadvantages of the feedback forcing approach.As discussed by (Iaccarino and Verzicco, 2003) the coefficients should be big enough for the restoring force to react on any frequency in the flow and thereby ensuring that the flow is independent of the coefficients.However, the governing equations tend to become very stiff when large values of K P and K I are used and therefore the method normally tend to suffer from being numerically unstable.Here we use K P = K p S FC where k p = 10 kgs − 1 m − 2 and K I = K i S FC where k p = 2 kgs − 1 m − 2 , which we found was close to the biggest values we could use and still obtain converged forces in all simulated cases.Besides proper selection of K P and K I , we also to some extent overcome the stability issue by solving the governing equations implicitly and using low relaxation factors (we use a relaxation factor of 0.4 for the momentum equations and 0.1 for the pressure correction equation).In addition, we find that numerical stability can be improved by slowly ramping up the IB forces in the beginning of the solution.Therefore, we usually multiply the IB force with a factor which grows linearly from 0 to 1 within the first 100 iterations.
The forcing parallel to the IB is due to friction and is computed by assuming that the near wall flow can be expressed locally via the logarithmic law of the wall for rough surfaces.Thus, we assume that the near wall velocity parallel to the wall, u p , is given by: Here d n is the normal distance from the velocity location to the IB, u * is the friction velocity, n o is the roughness height and κ = 0.4 is the von Kármán constant.
To utilize Eq. ( 2) for computing the friction force we extract the velocity parallel to the IB at the PP, u PP p using linear interpolation from the neighbouring cells and insert this velocity into Eq.( 2) with d n = d PP to obtain the friction velocity The friction force acting at the FP is then where ρ is the density and S FC is the area of the part of the IB which is inside the FC.Finally, the force F FP p is transferred to the computational grid in the same way as the normal force.

Boundary conditions for turbulence
In the present work we use the k − ω SST model as a closure for the RANS equations.By assuming local equilibrium between production and dissipation it can be shown that the following conditions apply for k (turbulent kinetic energy) and ω (specific rate of dissipation) in the logarithmic layer: where β = 0.09 is a constant and n is the distance to the FP.
These equations are used directly to set the value of k and ω, respectively in the centre of each FC.In practice we use a feedback approach similar to Eq. ( 1) to drive k and ω towards the desired values given by Eqs. ( 5)-( 6).

Computational domain and boundary conditions
The simulations are performed in a Cartesian domain as sketched in Fig. 5.
The domain has dimensions (L x , L y , L z ), where x is in the flow direction, z is vertical upwards and y is perpendicular to x and z.The number of cells in each direction is denoted (N x , N y , N z ).A near tree region of dimensions (l x , l y , l z ) is defined where the computational cells are cubic with side length Δ. Outside of the near tree region the cells are stretched towards the outer boundaries.
The boundary conditions are as follows: velocity prescribed according to the measurements at the upstream and top boundaries, zero gradient at the outlet, periodic conditions at the lateral boundaries and symmetry conditions (slip wall) at the bottom boundary.The reason for using symmetry instead of no-slip conditions at the bottom boundary is to avoid that the lower part of the imposed inflow velocity profile develops too much through the computational domain.The parameters governing the grids used for the RANS and DES are shown in Table 1.
The RANS simulations aim only at predicting the forces on the tree and therefore we use a high concentration of cells in the near-tree-region and allow a quite aggressive stretching towards the outer boundaries.On the other hand, the DES is focused on simulating the wake with a resolved turbulent inflow and therefore a milder stretching is used both up and downstream of the near-tree-region for these simulations.To limit the computational costs for the DES, the grid density is reduced in the near-tree-region compared to the RANS.The sensitive of the simulations to grid resolution is shown in the Appendix.

Inflow turbulence
The turbulent inflow used for the DES was generated using the spectral tensor model by (Mann, 1994) and introduced in a cross-section located 1.375H upstream of the tree using the technique proposed by (Troldborg et al., 2014).The advantage of this method is that it allows prescribing both the wind shear and turbulence levels according to a desired level but a drawback is that the flow is not fully in balance and hence undergoes some development as it moves downstream.
In order to resemble the measured conditions, we fit the parameters governing the spectral tensor in the algorithm by Mann to the measured spectra.A comparison of the measured and fitted spectrum is shown in Fig. 6.The used turbulence box has dimensions L X × L Y × L Z = 160H × 5H × 5H and cubic cells with side length H/25.6.
The fit is good for the small and medium turbulence scales but deviates at low wave numbers.The deviation is partly because the flow in the tunnel is not very representative of a typical neutral atmospheric boundary layer.It should be emphasized that establishing an exact match of the inflow conditions is further complicated in this case since the turbulence as described above is not fully in balance neither in the simulations nor the measurements and hence undergoes some development from upstream to downstream.In practice we overcome this issue by scaling the turbulence level such that the free-stream conditions at the wake positions fits the measured conditions at the same positions.

Results
As mentioned in Section 2, the 3D printed version of the tree has significant surface imperfections which are not present in the modelled geometry.These imperfections can only be modelled as an increased roughness.Therefore, we have conducted a study of the influence of tree surface roughness height on the predicted forces and the wake.In the following we will show the impact of n 0 on the forces whereas its impact on the wake is shown in the Appendix.Fig. 7 shows the simulated drag forces split into form and friction drag contributions at a fan speed of    100% for different surface roughness levels.As seen, the predictions are very sensitive to n 0 (especially for n 0 > 10 − 4 m) and show that it is important to accurately characterize the surface roughness to achieve accurate force predictions.In the following we will assume a roughness height of n 0 = 3⋅10 − 4 m because this was found to give a good agreement with the measured force at a fan speed of 100%.From a validation point of view, a smoother model tree would have been preferable since Fig. 7 suggests that the friction drag starts to play a minor role for roughness heights below 10 − 4 m (and especially when n 0 < 10 − 5 m).Nevertheless, given that real trees also have rough surfaces, the results in Fig. 7 point to that the assumption of negligible contribution by frictional forces used by (Endalew et al., 2009) and (Chester et al., 2007) may not always be valid.
It should be emphasized that the assumption of a uniform roughness distribution probably is not very adequate for the present tree but we have assumed it anyway since no estimates of the roughness distribution are available.
Fig. 8 compares the measured and simulated forces and corresponding drag coefficient as a function of u 0 , where u 0 denotes the freestream velocity at z/H = 0.74 (approximately mid crown height).The drag coefficient is computed based on the tree's stream-wise crosssectional area, A, and u 0 , i.e.: The agreement is good and suggests that when the surface roughness is chosen adequately then the present IB model can predict the forces on the tree rather accurately at different wind conditions.The drag coefficient is observed to be rather constant with wind speed and has a value which is in the same order as reported in earlier work (Rudnicki et al., 2004).The constant drag is expected because the tree is stiff and therefore is not streamlining with increased wind speed.Having shown that the model is capable of predicting the forces on the tree, the next step was to test its capability for predicting the wake flow.For this we use DES.
Initially we provide a qualitative impression of the simulated flow field by showing a snapshot of the simulated velocity and absolute vorticity contours in a vertical cross section at y = 0, i.e. through the stem of the tree.The plot is shown in Fig. 9.
The tree produces a deep and very turbulent wake which is characterized by small scale turbulence.The wake is deepest in the stem region but clearly expands upwards with downstream distance.The vorticity contours reveal the shear layers forming behind the individual branches and in some cases what appears to be vortex shedding although it is a bit difficult to see because of the strong ambient turbulence, which rapidly distorts the shed vortices.It should be noted that the simulations do predict vortex shedding from the stem region but it is not clearly visible in the present y = 0 plane.
In order to compare quantitatively the measured and simulated wake flow we gathered 3 second statistics.Since the measurements were carried out for T = 120 s this amounts to N 3s = 40 ensembles of 3 second long time series.The corresponding number of ensembles extracted from the simulations was only N 3s = 6.These simulated ensembles were extracted from two independent simulations, which had the same inflow statistics but used different turbulence seeds in the spectral tensor model by (Mann, 1994).Each simulation were run for a total time of 10.25 s but only the last 9 s were used in the analysis to avoid the startup transients and ensure a fully developed flow.Figs. 10 and 11 show the contours of the mean and standard deviation of the streamwise velocity at x/H = 0.63 and x/H = 1.90.Note that the velocity in both cases are scaled with the measured free-stream velocity u 0 at z/H = 0.74.In the near wake there is a good qualitative agreement between measurements and simulations.The wake deficit has similar shape with a slight asymmetry and the standard deviation is increased in the same regions which is most apparent at and above the upper part of the tree.In the far wake the comparison is less favorable.However, at this position, the wake is also less pronounced due to the high ambient turbulence and therefore it is harder to identify the perturbations induced by the tree.
Figs. 12 and 13 show the comparison of the measured and simulated mean and standard deviation of the streamwise velocity in the y = 0 plane through the tree.
The agreement between the measured and simulated mean velocity is good although the simulations consistently predict a larger deficit than observed in the measurements.However, the difference is almost always within the uncertainty in the measurements.The comparison of the standard deviations reveal a less favorable agreement and there is naturally more variations in the predictions.Whereas, the agreement at x/H = 1.90 generally is good at all heights, the predictions at x/H = 0.63 is noticeably lower than in the measurements especially at the lower positions.However, the simulation and measurements follow the same trend at x/H = 0.63 with a minimum in the standard deviation at around z/H = 0.5 and a significant increase in turbulence levels for z/ H > 0.75 with a peak at about tree height (z/H = 1).
An interesting outcome of the comparison is that both the predicted and measured turbulence level at x/H = 1.90 is lower than in the undisturbed flow.This is surprising but can be explained as an increased level of dissipation induced by the tree as discussed by (Ayotte et al.,   1999).It should be noted that the same phenomenon was observed in the full-scale experiment by (Dellwik et al., 2019).Although a one-to-one comparison cannot be made to their work because of differences in flow conditions and tree representation, it still suggests that the reduced wake turbulence is not an artifact of the flow imbalance in the WindEEE Dome.
Fig. 13 reveals that most of the simulated wake recovery occurs before x/H = 0.63 and that the turbulence level reaches a plateau beyond this point.This could indicate that the above mentioned deviations in turbulence between measurements and simulations at x /H = 0.63 is due to a faster initial turbulence decay in the simulations.Note that the small kink observed in both mean and standard deviation at approximately x/H = 0.3, z/H = 0.61 is caused by the presence of a branch close by.
To give more insight to the behavior of the turbulence levels, Fig. 14 presents the spectral characteristics of the streamwise velocity at z/H = 0.35 and z/H = 1.0.Note that the reason that the simulated spectra are noisier than those from the measurements is that they are based on fewer ensembles of time series as mentioned above.
From Fig. 14(a) it is clear that the increased standard deviation in the measurements at (x/H, z/H) = (0.63, 0.35) is caused by an increase of energy at frequencies above 10 Hz.Fig. 14(b) reveals that the simulations predict a similar increase at these frequencies and that the lower standard deviation observed in the simulations at this position is caused by a reduction in the energy for frequencies below 1Hz.It is unclear what is causing this decrease in energy in the simulations but it cannot  be explained from a non-balanced turbulence field because the decrease is not observed when the tree is not present.
Apart from this deviation, the agreement between the measured and simulated spectra is generally rather good at all locations.

Discussion
The above results show that with a suitable surface roughness, the IBM is capable of predicting the aerodynamic forces on a tree skeleton rather accurately.A shortcoming of the present validation is however that we do not have a proper characterization of the tree's surface roughness and therefore we do not know whether the roughness height used in the modelling is actually representative of the tree's true surface roughness.Based on a visual inspection of the tree surface we did however expect a roughness height in the order of 10 − 4 m prior to performing the simulations.Given that real stem and branch surfaces often are rough either due to the structure of the bark or to imperfections in the tree geometry, simulations using very low roughness values or noslip conditions may underestimate the drag force.Another drawback of the present study is that the non-balanced flow in the wind tunnel makes the validation of the numerical model uncertain.Although the validation is not perfect, the above results still clearly shows that the model can capture the main physics governing the forces on the tree and the flow in it's wake.The sensitivity to tree surface roughness was higher than expected and shows that at least for some types of trees it is important to accurately characterize the roughness to achieve accurate force predictions.However, in practice it is very difficult to estimate the surface roughness for complex structures like trees as well as to make a consistent split between the parts of the tree that is resolved and the parts that need to be modelled as increased roughness.Therefore it is also a parameter which is rarely known for trees and hence challenges the use of the IBM for predicting tree forces in general.Nonetheless, the dependence on roughness is still a major advancement compared to models which requires specification of a drag coefficient because the roughness is a property of the surface and not the whole tree as in (Dellwik et al., 2019).Even though simulations of trees with leaves inevitably will require specification of leaf area density and drag coefficients, the combination with the IBM for simulating the stem and main branches should give a better prediction of the wake especially in the stem region.
The method presented here could be expanded to include the unresolved branches, twigs and leaves in a similar manner as the IM method suggested by (Endalew et al., 2009).However, the immersed boundary method presented here has the advantage over their approach with a surface-resolving grid, since it more easily can be applied for different complex geometries.
The method of using the spectral tensor model by (Mann, 1994) to generate the inflow turbulence for the simulations is essential for the present study because the spectral tensor can be fitted to the actually measured spectral characteristics and thereby allows mimicking the non-ideal inflow conditions in the WindEEE Dome.Since in practice ideal inflow conditions hardly ever exist, we argue that our simplified representation of the ambient flow field is a promising way forward.

Conclusion
Measurements and simulations of the forces on a leafless tree and the wind field in it's wake have been presented.In the simulations the tree is represented by the immersed boundary method (IBM), while the flow field is represented with either Reynolds Averaged Navier Stokes (RANS) or detached eddy simulations (DES).The agreement between measurements and simulations was in general found to be good both in terms of forces and wake flow.However, the good agreement in forces was obtained only after calibrating the surface roughness applied in the simulations to the measured force at a given wind speed.This calibration was necessary because the "true" roughness characteristics of the tree was unknown.The simulations, however revealed a strong dependence of the forces to the modelled surface roughness of the tree and therefore also suggest that the roughness is an important parameter to quantify in the characterization of trees.

Declaration of Competing Interest
None.
plots.As expected an increase in surface roughness leads to a deeper deficit and increased turbulence levels in the wake.The wake characteristics are only slightly affected by a change in roughness for n 0 < 10 − 4 m, which is consistent with the corresponding sensitivity of the drag force (see Fig. 7).

Sensitivity to Reynolds number
The model tree tested in the wind tunnel is at a scale of approximately 1:16 compared to the real tree, which resulted in a rather low Reynolds number based on tree height of Re H = 1.75⋅10 5 .To investigate any potential issues relaed to the down-scaling, we performed DES at a Reynolds number which is 16 times larger than in the simulations presented in Section 5. Fig. A3 compares the mean and standard deviation of the streamwise velocity obtained at the two different Reynolds numbers, while Table A1 compares the same statistics for the drag.The uncertainty of the drag statistics is expressed in terms of the standard deviation of the mean.As seen the results obtained at the different Reynolds numbers are very similar and therefore we conclude that the results obtained in this paper are not significantly affected by scale effects and that our IBM therefore should also be valid on full scale.

Table A1
Mean and standard deviation of drag force at different Reynolds numbers.The standard deviation of the mean (σ/ ̅̅̅̅̅̅̅ N 3s √ ) is also indicated .

Fig. 1 .
Fig. 1. 3D printed tree in the WindEEE Dome (left) and a close up of the tree skeleton (right).

Fig. 2 .
Fig. 2. Sketch showing (using black dots) the locations where wake measurements were recorded; a) Side view; b) Top view.

Fig. 7 .
Fig. 7. Simulated forces as a function of tree surface roughness.The used roughness height is indicated with a dashed line.

Fig. 8 .Fig. 9 .
Fig. 8.Comparison of measured and simulated drag force (left) and drag coefficient (right) as a function of free-stream velocity at z/H = 0.74.

Fig. 10 .
Fig. 10.Contours of streamwise velocity at x/H = 0.63 (top) and x/H = 1.90 (bottom) downstream of the tree.Left plots are measured and right plots are simulated.

Fig. 11 .
Fig. 11.Contours of standard deviation of the streamwise velocity at x/H = 0.63 (top) and x/H = 1.90 (bottom) downstream of the tree.Left plots are measured and right plots are simulated.

Fig. 12 .
Fig. 12. Vertical profiles of mean (left) and standard deviation (right) of the streamwise velocity in the wake of the tree.Simulations are shown as full lines and measurements as dots.The standard error of the estimates are indicated as shaded areas for the simulations and as error bars in the measurements.

Fig. 13 .
Fig. 13.Mean (left) and standard deviation (right) of the streamwise velocity as a function of x at different heights.Simulations are shown as full lines and measurements as symbols.The standard error of the estimates are indicated as shaded areas for the simulations and as error bars in the measurements.Note that due to the scale on the y-axis the error bars in the measurements are barely visible.

Fig. A1 .
Fig. A1.Predicted streamwise force as a function of number of cells per tree height.

Fig. A2 .
Fig. A2.Mean and standard deviation of the streamwise velocity at x/H = 0.63 (top) and x/H = 1.90 (bottom), respectively for different levels of tree surface roughness.

Fig. A3 .
Fig. A3.Mean and standard deviation of the streamwise velocity at x/H = 0.63 (top) and x/H = 1.90 (bottom), respectively for two different Reynolds numbers.

Table 1
Parameters defining the meshes used for the RANS and DES.H is the tree height.