Stems were approximated as truncated cone meshes, represented using quadratic tetrahedral elements. Each (fully grown) mesh was divided into elements covering one fifth the radius, one fifth the height and one sixteenth of the circumference. Data presented in \citet{waghorn_influence_2007} was used to create a regression to estimate stem radius (\(r_{s}\)) from stocking rates based on the slenderness ratios. The regression presented in Equation \ref{eq:slenderness} was calculated from reported slenderness ratios for stocking rates of between 275 and 1457 stems per hectare and assuming a stem height of 15 m.

Where \(S\) is the stocking rate, assumed to be either open grown, one stem per hectare or 741 stems per hectare.

As trees grow, the newest cells are superimposed on the existing structure, creating a time evolving stress state. In order to accommodate these time dependent changes, the mesh needs to be ’grown’, here time is treated in discrete steps. Wind effects are not accounted for during growth, however, self weight and growth stress are. Growth stresses develop as cells divide in the cambium, resulting in the new cells being in a non-stressed state on the surface of a pre-stressed structure. During secondary wall formation there is a tendency for cells to contract longitudinally, for more details see \citep{archer_growth_1986}. As the pre-existing structure restrains further contraction a stress profile of longitudinal compression at the centre of the stem and tension at the outer edge forms.

The pre-stressed state of the cambium from gravity and growth stresses produces the need to successively ’grow’ the mesh which is to be used so that each mesh addition is added in a non stressed state on top of a stressed surface. Mesh growth is achieved by first defining an initial state, this state can be thought of as a seedling. The seedling is subjected to gravitational forces, deformation occurs and an updated mesh is created. The positions of the nodes and vertices of the deformed mesh are used to calculate the positions of new vertices and nodes to be added in order to represent cambial growth in a non stressed state onto the stressed surface along with apical growth to a predefined height. This procedure is repeated for five time steps imposing gravity (and in required cases growth stress), regardless of the final age, i.e. a tree with wind imposed at 5 years utilises 1 year steps, a tree with wind imposed at 15 years has 3 year time steps.

The material constants for Pinus radiata reported by \citet{Davies_2016} were used to parametrise the model. Four sets of constants representing extreme values of green density and stiffness for the species were used. Note the abbreviations of the sample type used here, the first letter defines the density as either \(H\) (High density) or \(L\) (Low density) and the second letter defines the MFA as either \(W\) (Wide) or \(S\) (shallow).

Linear interpolation was used to approximate material constants between two samples to provide a gradient of material properties from the corewood to the outerwood. The TRP of cellular properties was investigated using the gradient of elastic constants and Poisson ratios from different samples by applying simulated wind loads to stems with different radial profiles.

While the TRP of Pinus radiata usually follows a non-linear increase for density and a non-linear decrease for MFA \citep{burdon_juvenile_2004}, during the first 15 years the trend can be approximated as linear, after this however the non-linear behaviour becomes more prevalent as the properties stabilise. It should be noted that there is much variation between individuals and this is not universal. The TRP was investigated by calculating the appropriate elastic and shear moduli along with Poisson ratios for each point in the stem. The linear interpolation was calculated between two samples chosen to represent corewood and outerwood, Equation \ref{eq:constants} shows the general form.

Where \(r\) is the stem radius at the current point, \(r_{t}\) is the total radius at 15 years, \(λ_{i}\) and \(λ_{o}\) are the values of the material properties obtained from \citet{Davies_2016} at \(r=0\) and \(r=r_{t}\) respectively.

Due to the micro-structure of wood the native coordinate system for describing material properties is not the same as the global coordinates system used to impose external forces on the stem. Transformation between the two systems is needed. The \(rtl\) local system used for experimental work presented in \citet{Davies_2016} is used in combination with the interpolations presented in Equation \ref{eq:constants} and presented as the stiffness matrix in Equation \ref{eq:stiffness_mat} using Voigt (engineering) notation at any point in the stem. Because the stiffness matrix is calculated in the local coordinates it needs to be converted into an \(xyz\) system in order to apply wind loadings in a sensible fashion, Equation \ref{eq:trans_mat} provides the transformation matrix. The use of Voigt notation allows for transformation of the elastic constants into the global system, via Equations \ref{eq:cm_trans} to \ref{eq:stress_tans}. It is assumed that there is no spiral grain occurrence within the stem, the local \(l\) axis is always parallel with the global \(z\) axis. Material constants in the longitudinal direction are parallel with the \(z\) axis, i.e. there is no correction angle applied to account for taper.