Finite element analysis of trees in the wind based on terrestrial laser scanning data

Wind damage is an important driver of forest structure and dynamics, but it is poorly understood in natural broadleaf forests. This paper presents a new approach in the study of wind damage: combining terrestrial laser scanning (TLS) data and finite element analysis. Recent advances in tree reconstruction from TLS data allowed us to accurately represent the 3D geometry of a tree in a mechanical simulation, without the need for arduous manual mapping or simplifying assumptions about tree shape. We used this simulation to predict the mechanical strains produced on the trunks of 21 trees in Wytham Woods, UK, and validated it using strain data measured on these same trees. For a subset of five trees near the anemometer, the model predicted a five-minute time-series of strain with a mean cross-correlation coefficient of 0.71, when forced by the locally measured wind speed data. Additionally, the maximum strain associated with a 5ms −1 or 15ms -1 wind speed was well predicted by the model (N=17, R 2 =0.81 and R 2 =0.79, respectively). We also predicted the critical wind speed at which the trees will break from both the field data and models and find a good overall agreement (N=17, R 2 =0.40). Finally, the model predicted the correct trend in the fundamental frequencies of the trees (N=20, R 2 =0.38) although there was a systematic underprediction, possibly due to the simplified treatment of material properties in the model. The current approach relies on local wind data, so must be combined with wind flow modelling to be applicable at the landscape-scale or over complex terrain. This approach is applicable at the plot level and could also be applied to open-grown trees, such as in cities or parks.


Introduction
Wind damage to forests is a significant but poorly understood driver of the terrestrial carbon cycle (Espírito-Santo et al., 2014). Trees respond to strong winds by leaning and swaying in a dynamic manner. This response has been characterized using static models (Hale et al., 2015), dynamic models with simplistic tree geometry (Kerzenmacher and Gardiner, 1998) and dynamic models which represent the effect of branches . The literature has mostly focused on wind damage in conifer plantations due to their importance to the forestry industry (Dupont et al., 2015;Hale et al., 2015Hale et al., , 2012. In broadleaf forests the variety of tree sizes and shapes influence their response to wind forcing, meaning that finite element analysis is necessary to model the detailed effects of wind on trees Rodriguez et al., 2008;Spatz and Theckes, 2013). Previous finite element analysis studies have been limited to small numbers of trees due to the difficulty of manually mapping tree geometry (Moore and Maguire, 2008;Sellier et al., 2006). However, terrestrial laser scanning (TLS) data offers the potential for dramatic advances in this area. Recent developments in the extraction of accurate 3D tree models from TLS data can provide detailed geometrical information for a scanned tree in a matter of minutes (Raumonen et al., 2015;Calders et al., 2015;Gonzalez de et al., 2018).
We used 3D tree models derived from TLS data as the basis of a finite element analysis, applied using the engineering software Abaqus (Simulia Software Company, 2017). This approach allowed us to simulate the dynamics of broadleaf trees without the need to manually map their geometry. In addition to its application in natural forests, this approach could be useful in testing proposed interventions on urban amenity trees, or valuable veteran trees, in a model environment before implementation (Reiland et al., 2015).
This paper tests the hypothesis that finite element models developed using TLS data can be successfully used to simulate the behaviour of broadleaved trees under applied wind loading. We describe the method and discuss the assumptions involved. We also provide software to implement this technique. We test the model against field data collected in Wytham Woods, UK, for which we have both TLS data and field measurements of strain in the tree stem and wind speed. This validation takes three forms: 1 We predict the exact time-series of strain produced during a fiveminute period based on high time-resolution wind data. 2 We calculate the maximum hourly strain produced on a tree given maximum hourly wind speeds of 5 ms −1 or 15 ms -1 . We also estimate the critical wind speeds (CWS) at which trees will break by extrapolating from long term field data (Hale et al., 2012). 3 We predict the fundamental frequencies of the trees and compare these to the frequencies extracted from the field data using Fourier analysis.
Finally, we discuss the results of this study in the context of other work in the field.

Materials and methods
The workflow involved field data collection and analysis, TLS data collection and modelling work, followed by three separate tests of the model against the field data ( Fig. 1).

Study site and field data
Wytham Woods (51°46′27.2″ N 1°20′20.1″ W) is a mature temperate deciduous woodland near Oxford, UK. The woodland contains approximately 950 trees per hectare and the dominant species are sycamore (Acer pseudoplatanus), ash (Fraxinus excelsior), beech (Fagus sylvatica) and English oak (Quercus robur). We measured the bending strains produced on 21 trees (eight ash, nine sycamore, four birch (Betula spp.) and one oak) ranging from 12 to 24 m tall. One of the sycamore trees had four main stems below 1.3 m, each of which are treated as separate trees in the forest census and were instrumented as such in this study, we therefore refer to them as separate trees in this paper. Seven of the trees were located near a canopy walkway structure, and the remaining 14 trees were located 230 m away. Average canopy height was 21-23 m in both locations.
Strain was measured using the method developed by Moore et al., (2005): pairs of strain transducers were fixed to the trees at approximately 1.3 m and perpendicular to each other. These transducers were then connected to Campbell Scientific (Logan, UT, USA) data loggers and measured at 4 Hz (two CR23xs and three CR1000s). The limited storage capacity of the CR23x data loggers was supplemented using a Raspberry Pi (Cambridge, UK).
In the field, strain gauge resistances drift with environmental factors. We partially corrected these in situ using a variable resistor within a wheat-stone bridge circuit. We further accounted for this drift by subtracting a 5000 point (20 min) running mode from the signal. This is based on the assumption that the tree will regularly pass through its equilibrium position. This same technique was used in Wellpott (2008) and Hale et al. (2012Hale et al. ( , 2015 and we found our results to be relatively insensitive to this window length (SI-1).
Wind speed was measured in winter using three cup anemometers (A100LK/5 M, Vector Instruments, Rhyl, Wales) mounted at 5 m, 10 m and 15 m heights above ground level on a canopy walkway next to the focal trees and supplemented, notably during the 5-minute high-resolution validation period, with a 3D sonic anemometer (Campbell Scientific CSAT3) mounted at 15 m. In summer, we used a Gill (Gosport, UK) Sonic-1 anemometer recording at 1 Hz mounted at 15 m. Local climate data are available through a meteorological station approximately 1 km away, just outside the forest, operated by the Environmental Change Network (ECN, www.ecn.ac.uk). This station records mean hourly wind speed as well as the maximum 5 s mean within that hour (5 s maximum gust speed). The wind data taken from anemometers within the canopy were considered unrepresentative, especially in summer, because of localized sheltering effects. We therefore used the ECN meteorological station's hourly maximum wind speeds for all critical wind speed analysis and select the maximum strain on each tree within each hour. A comparison of the different measures of wind speed is given in SI-2.

Development of quantitative structure models
TLS data were collected in Wytham Woods, Oxford, as part of a separate project (Calders et al., , 2016. Trees were manually extracted from the point cloud and cylinder fitting algorithms, developed primarily for biomass estimation, were used to extract the tree shapes (Åkerblom, 2017;. This method represents each tree as a series of cylinders -a quantitative structure model (QSM).
We adapted these QSMs by removing cylinders under 2 cm diameter, which are generally the least accurate in the QSM reconstruction due to limitations of TLS scanning resolution (see Fig. 2). We then longitudinally averaged neighbouring cylinders, replacing each pair of neighbouring cylinders with a single cylinder with the average length, radius and direction of the original pair (Åkerblom, 2017). This reduced the number of cylinders from approximately 4000 per tree to 500, without significantly affecting the total woody volume of the tree. This process resulted in a new set of QSMs comprised of cylinders with a higher length to radius ratio, allowing us to use beam theory to build a dynamic model of the tree without significantly altering the overall structure (Carrera et al., 2011). It is important to note that the Wytham Woods TLS data were collected in ideal conditions based on a 10 m grid using a Riegl (Horn, Austria) VZ-400 with 0.04°angular resolution (Wilkes et al., 2017). Most of the trees were scanned in winter, but those near the walkway were scanned in summer and the leaves removed through filtering (Vicari, 2017). Optimization of the QSM approach shows that there is very little variance between the QSMs in  these conditions .

Development of finite element models
To build a finite element model in Abaqus, each cylinder in the QSM was represented as a two node Euler-Bernoulli beam element (B31) and the orthotropic features of wood were neglected (Sellier et al., 2006). Since we intend this approach to be applicable wherever TLS data is available, we did not measure the material properties of the trees individually. Values of green wood density, , gw elasticity E , , gw and modulus of rupture, MOR gw , were taken from the literature (Lavers, 1983;Niklas and Spatz, 2010) (Table 1). Breaking strain was defined as = MOR E / break gw gw (Niklas and Spatz, 2010). The intra-specific variation in elasticity is between 14 and 20% (defined as standard deviation divided by mean value) in sycamore, ash and oak. No estimates of variation in green wood density were given. The inter-specific variation in wood density is 20% and elasticity 25% in Wytham Woods (defined as the range divided by the minimum value, Table 1). The values used in finite element analysis are given in Table 1.
Root anchorage was modelled as a universal joint with two rotational degrees of freedom and symmetric elasticity = E root 5 MPa in each direction. Energy dissipation at the root-soil boundary was modelled by a dashpot with coefficient = C 10 s . Previous authors used tree pulling data to tune these model parameters, which gives accurate results in the case of a few trees (Sellier et al., 2006;Sellier and Fourcaud, 2009). However, we are attempting to generalize to many trees without intensive field measurements, so the root parameter values were chosen because the model is insensitive to these parameters in this range (Table 2). Roots are therefore included in the model primarily so that future studies can tune these parameters if the necessary data are available.
Applying gravity in a realistic way proved difficult, since the trees were scanned in their gravity deformed positions. Branches develop reaction wood and grow in irregular shapes to counteract the force of gravity (Scurfield, 1973), whereas our model assumes cylindrical branches with isotropic material properties. We tested three approaches to modelling the effect of gravity and found a small but interesting difference between them (see . For all analysis in the main text of this paper we first applied an inverse gravity force to the model tree, extracted the nodal displacements and then restarted the analysis applying a realistic gravitational force to the displaced geometry. This is the most physically realistic approach since branches are not displaced before the analysis begins and the additional load due to self-weight of the displaced crown is accounted for.

Dynamic simulation
The dynamic simulation of a tree's response to wind forcing was conducted using the Abaqus Aqua module (Simulia Software Company, 2017). Each tree was modelled individually, meaning that the dynamic response does not include crown clashing which is thought to be significant in some forest structures (Rudnicki et al., 2003(Rudnicki et al., , 2001. In all simulations, we used a non-linear solution method to account for large displacements in the tree structure. The equation of motion takes the form where q q , and q are the nodal displacement, velocity and acceleration vectors. M, D and K are the mass, damping and stiffness matrices respectively, and G and F D represent the gravitational and wind induced forces. The output variable, strain, is calculated internally by Abaqus and we use the strain at the approximate position of the strain gauges (1.3 m), so as to correspond with the field data.
While horizontal wind velocities were fully represented, the vertical profile of wind speed was limited to a power law fit (Eq. 2), using field data from three heights within the canopy ( Figure SI-2).
where v is the wind speed at a given height h, v 0 is the reference wind speed at reference height h ref and is a variable parameter determined by the local wind profile ( Figure SI-2). The wind drag force, F D , per unit length on each cylinder of the QSM is given by where Amp is a user defined amplitude that describes the horizontal variation in wind speed, air is the density of air, C D is the drag factor, r is the radius of the beam, and V fn is velocity difference between the wind and the beam. As in previous work (Ciftci et al., 2013;Moore and Maguire, 2008;Sellier et al., 2006;Sellier and Fourcaud, 2009), we assumed a constant drag factor, although in fact the drag factor varies with wind speed and with the size of the branch. The most significant source of damping, that due to aerodynamic drag, was modelled implicitly in the Aqua module. Material damping was modelled using Rayleigh's hypothesis (Rayleigh, 1877).
where M W is the mass proportional damping and K is the stiffness proportional damping. Following Sellier et al. (2006) we assumed stiffness proportional damping would be insignificant in the high Fig. 2. Left -QSM directly from the fitting process. Right -simplified QSM with longer cylinders and fewer small branches -this is our input for finite element analysis.

Table 1
Material properties of green wood used in finite element analysis. Values in brackets are standard deviation and number of samples tested. Quercus rubra was used as a substitute for Quercus robur, since data for the latter was not available from the same testing method. This affects only a single tree in the frequency analysis validation.
Reynolds number regime in question, and so set = 1 × 10 −3 . Contact between branches was allowed but no attempt was made to accurately model the energy transfer of these contacts.
Multiple simulations were carried out to test the sensitivity of the model. All continuous parameters were tested at ± 20%. Non-continuous parameters, such as the simplification of the QSM, were tested at different levels.

Validation techniques
All the strain data used for validation were collected in winter, when the absence of leaves makes the modelling and interpretation simpler. Future work will parameterize the effects of leaves, but the change in Reynolds number and drag coefficient make the full treatment of wind flow over leaves highly complex (de Langre, 2008).

Time series validation
For the time-series validation, we use a five-minute subset of data from 14:00-14:05 on 28 th March 2016, chosen for high wind speeds and high (10 Hz) wind measurement resolution. Two trees were removed due to problems with the field data, in the first case the balancing circuit failed and in the second the physical connection to the tree came undone, leaving five trees with reliable data. We therefore compared measured and modelled strain time-series for five trees close to the anemometers. We did not expect an exact correspondence for the 5minute time-series, since the wind forcing applied to each model tree was identical and reflects the wind regime measured at the position of the canopy walkway. In the field, each individual tree experiences a slightly different wind forcing due to its position (Kerzenmacher and Gardiner, 1998). A previous finite element analysis study found a reduced correlation at 12 m from the wind measurements, although no correlation statistics were given (Sellier et al., 2008).

Hourly maximum strain validation
In this validation, we used the five trees from time series validation as well as twelve additional trees situated approximately 230 m away. Two of these additional trees were excluded due to a poor connection between the tree and the strain gauge. We used hourly maximum wind speed measurements taken from the ECN meteorological station. For each hour we extracted the maximum strain value for every tree. All but four trees had over 500 h of data available, and for these four trees we had 165 h of data. We tested a robust maximum selection technique, using the mode of the maxima of subsections of the hours data, but found no increase in coefficient of determination nor reduction in root mean squared error . Details of the strain data processing and reference to the software and data repository to repeat this analysis are given in SI-1. We regressed the hourly maximum strain data against the squared hourly maximum wind speed and so predicted the strain at 5 and 15 ms −1 wind speeds. We used this regression line to extrapolate the relationship up to the breaking strain, giving an estimate of the critical wind speed from the field data. In the simulation, we artificially increased the wind forcing until the lowest beam in the trunk reached its breaking strain.

Fundamental frequency validation
We used a Welch's power spectral density function (Krauss et al., 1994) to extract the fundamental frequency for each hour of field data during winter and take the mean. This analysis only requires one strain gauge per tree, so three of the four strain gauge failures had no effect here, meaning that we can carry out this analysis for 20 trees. Using Abaqus, we predicted all the undamped sway modes from the QSMs using a subspace method (Viberg, 1995). This can be performed separately from the dynamic analysis and is much faster, it is also available in the free version of Abaqus.

Time-series validation
The modelled five-minute time-series of strain was similar to that Fig. 3. Above -TLS data showing trees near the canopy walkway. The coloured trees are the ones used in this analysis. The black arrow shows the prevailing wind direction and is located at the point of the wind measurement. Left -simulated strain (coloured lines) against field data (black lines, inverted for clarity). The time-lag (field strain -simulated strain) as well as the cross correlation co-efficient are given for each signal. measured in the field (Fig. 3, mean cross-correlation = 0.71). We found the lag between measured and simulated strain was partially explained by the position of the trees in relation to the wind measurements. However, all modelled trees displayed some peaks that were not evident in the field data, and vice-versa. This was presumably due to differences in wind at the anemometer position and at the position of the tree. This time-series validation is similar to that given in Moore and Maguire (2008) and Sellier et al. (2008). It is useful to demonstrate that the modelled response resembles that of the actual trees, but only applicable at short distances from the anemometers due to the spatial variability of wind speeds.

Hourly maximum strain validation
This section examines the model's accuracy in predicting the strains produced by wind speeds of 5 and 15 ms −1 , and in extrapolating up to the critical wind speed (Fig. 4). It is important to note that, in the case of CWS, both the x and y axes are estimates, since the CWS was estimated from field data by extrapolating beyond the available measurements. Therefore, this is not a straightforward validation but rather a demonstration that predictions from field and model approaches are correlated. We used hourly maximum wind speeds (highest 5 s mean per hour) measured at the nearby meteorological station, since this avoids the need for a multiplication factor (gust factor). A more detailed discussion of the data analysis choices relating to the wind data is given in SI-2.
The model predicted the strain produced at 5 ms −1 with an R 2 = 0.81 and concordance correlation coefficient (CCC) of 0.91. At 15 ms −1 wind speed the fit was slightly reduced R 2 = 0.79 and CCC = 0.89. CWS predictions from field data and the model displayed a worse agreement, R 2 = 0.40 and CCC = 0.66, presumably because the model behavior departed from a pure square law before reaching the CWS (due to streamlining), while the best fit line derived from the field data does not. We saw no departure from a square law in our field data, but we did not have measurements at high enough wind speeds to test this.

Fundamental frequency validation
We found that finite element analysis predicts the trend of fundamental frequency (f 0 ) correctly (R 2 = 0.38 and CCC = 0.65, N = 20). Using sample mean material properties, instead of the species-specific material properties which are often unavailable, only decreased the R 2 value by 0.05. Interestingly, there is a systematic under-prediction of f 0 across the 20 trees (Fig. 5). This is likely due either to our simplistic treatment of material properties, which are constant throughout the tree, or our arbitrary trimming of all branches with a diameter under 2 cm. As noted in the figure legend, the fit is radically improved by removing the two birch trees whose fundamental frequencies were over predicted.
The free version of Abaqus can perform this analysis for any tree model consisting of fewer than 1000 beams. As previously noted this method predicts a dense distribution of modes about the main frequencies (Rodriguez et al., 2012), representing the slightly different shapes or directions in which the tree can sway. In the field data, this takes the form of a wide resonance peak, instead of distinct modes.  Field measured fundamental frequency against modelled fundamental frequency for 20 trees (one tree was excluded due to errors in the simulation). The goodness of fit is highly sensitive to the two birch trees for which the fundamental frequency was overestimated. Removing these two trees increases the goodness of fit (R 2 to 0.66, CCC = 0.83).

Sensitivity analysis
In this section, we examine the sensitivity of the model to the input parameters. We used three output variables as indicators of the models' response to a prescribed change in input parameters, the strain at a wind speed of 15 ms −1 , the critical wind speed estimate and the fundamental frequency ( Table 2). The model was highly sensitive to the details of the QSM used as input. We found that the simplification step, used to reduce the number of cylinders in the original QSM for finite element analysis, has a large effect on the result. The second largest sensitivity was due to material properties, especially wood elasticity. Natural variation in wood elasticity, where data are available, is between 14% and 19%, which is comparable to the range of our sensitivity analysis. Unfortunately, green wood elasticity is not often measured and different measurement techniques are often not directly comparable (Ruel et al., 2010).
The beam type is a property of the finite element analysis and determines the solution method used to calculate the strains in each beam (Simulia Software Company, 2017). The choice of beam type is one of the smaller sources of sensitivity in the model. The drag factor, C D , is an experimentally determined parameter which determines how much force is transferred to an object for a given wind speed (see Eq. 3). We treated the drag factor as a constant, following (Ciftci et al., 2012;James, 2010;Moore and Maguire, 2008;Sellier and Fourcaud, 2009), whereas in fact it changes with the Reynolds number of the flow (de Langre, 2008). These changes are small in the low wind speed regime of our measurements and several projects are currently investigating the changes in the drag factor for trees at high wind speed. The model is insensitive to root parameters in this range, which are poorly constrained in the literature.

Discussion
Overall, finite element analysis using QSMs from TLS data is a useful method for modelling the behaviour of trees in the wind (Table 3). It can predict the time-series of strain produced in the trunk of a tree, if the model is forced by wind speed data measured nearby. It can also predict the maximum strain for a given maximum wind speed, the critical wind speed at which the tree will break and the fundamental frequency of a tree. Since TLS data are becoming increasingly available for trees from a wide range of environments from tropical forests to city parks, this technique could now be used to model the dynamic behavior of a wide range of trees. This will allow us to revisit questions about the relationship between tree architecture and dynamics with realistic tree shapes instead of theory-based fractal or modelled trees. In particular, we could explore the relative importance of variations in material properties and tree architecture (Sellier and Fourcaud, 2009) and the differences between open grown trees and those in forests Kane et al., 2014).
The characteristics of the high time-resolution response of a tree to wind are important, since it is at these scales that dynamical amplification and damping processes occur, which alter the likelihood of wind damage (Ciftci et al., 2013;James et al., 2006;Spatz et al., 2007). However, wind speeds are highly spatially variable, so an array of anemometers is needed to represent the wind forcing. The variability of wind speeds can broadly be described by coherent gusts of approximately one tree height in the cross-stream direction and significantly longer in the stream-wise direction (Finnigan, 2000). In the current model, the wind regime at the position of each tree must be provided and we assume these wind regimes are identical. In the high-time resolution validation this led to some gusts hitting the trees which were not measured by the anemometer and vice versa.
The static response of the QSMs to an applied wind forcing was similar to the hourly maxima measured in the field. This is in agreement with recent work which demonstrated that, for Pinus sylvestris, a large part of the tree response was due to simple static effects (Schindler and Mohr, 2018). Within the range of wind speeds measured in this study the correlation between model and field data exceeded R 2 = 0.79, CCC = 0.89. As the wind speed increased this correlation decreased, possibly because the model accounts for non-linear effects due to the displacement of the crown while the extrapolation from field data did not. The critical wind speeds were systematically overestimated by the simulation. It is important to note that there are likely to be further nonlinear effects at high-wind speeds that are not accounted for in either the field data or the model. These are associated with changes in the drag factor, streamlining of the tree crown and possible plastic deformation of the wood. However, we have no data with which to quantify these effects and assume they are small in comparison with the main effect of increasing wind speed.
The fundamental frequency is an important parameter since a lower fundamental frequency means that energy is more easily passed from the wind to the tree and the tree is at a higher risk of damage (Schütz et al., 2006).The model predicted the fundamental frequency of trees with R 2 = 0.38, CCC = 0.65, but displayed a clear negative offset. This offset is likely due to either our simplistic treatment of material properties or the simplification of the QSM. We do not possess independent validations of the accuracy of the QSM reconstructions and therefore cannot separate errors due the QSMs and those due to the treatment of  material properties. A review of the material properties of trees worldwide (Niklas and Spatz, 2010) found a clear relationship between wood density and elasticity at the species level, but the potential for wide variation within species, especially in those selected for rapid growth, and within individual trees remains. Our sensitivity analysis shows that wood elasticity plays a larger role than wood density in driving tree dynamics, but it is also harder to measure (Ruel et al., 2010). Previous authors increased the accuracy of the finite element analysis by adjusting the material properties in order to fit the model to static field tests performed on the tree (Sellier et al., 2006). This method could prove highly useful for small numbers of urban or valuable veteran trees where proposed treatments, such as removing a large branch or adding supports, can be tested in a model environment before being carried out. Another possible addition to the model would be a more accurate representation of leaves and twigs, which must influence the motion of the tree. This is currently limited by the difficulty of retrieving small twigs and leaves from TLS data, but a parameterization could possibly be developed using summer (leaf-on) data in combination with the winter (leaf-off) data.

Conclusion
Previously, finite element analysis of trees has been confined to a small number of relatively short or even idealized trees. We demonstrate that combining TLS data and finite element analysis allows us to model the detailed dynamics of realistic and complex broadleaf trees without manually mapping their 3D architecture. We show that this approach can model the strains produced in a tree for a given timeseries wind forcing. Finite element analysis can also predict the effects of a static wind forcing, which was found to follow a similar trend to the hourly maximum strain -wind speed relationship. It can also predict the correct trend of fundamental frequencies for a set of trees given their geometry, although there is a systematic offset.
Several uncertainties remain, in particular relating to the accuracy of the TLS data in representing 3D tree architecture, the mechanical properties of roots, and the change in wind-tree interaction at high wind speeds. This approach is applicable in natural forest plots and would also be valuable in predicting the dynamics of park or city trees, and how they will respond to a proposed intervention.

Data availability
All the field data collected for this study is available online (https:// doi.org/10.5285/533d87d3-48c1-4c6e-9f2f-fda273ab45bc)as are the data processing scripts (https://github.com/TobyDJackson/ WindAndTrees_Wytham). The matlab scripts that convert QSMs to Abaqus input files are also available online (https://doi.org/http:// doi.org/10.5281/zenodo.894543). We used the commercial finite element analysis software Abaqus (Simulia Software Company, 2017), but a student version can be freely obtained on which a more limited analysis can be run -this is clearly explained in the software readme files.