EAGLE model¶
This section of the documentation gives a brief description of the different components of the EAGLE subgrid model. We mostly focus on the parameters and values output in the snapshots.
Entropy floors¶
The gas particles in the EAGLE model are prevented from cooling below a
certain temperature. The temperature limit depends on the density of the
particles. Two floors are used in conjonction. Both are implemented as
polytropic “equations of states”\(P = P_c
\left(\rho/\rho_c\right)^\gamma\) (all done in physical coordinates), with
the constants derived from the user input given in terms of temperature and
Hydrogen number density. The code computing the entropy floor
is located in the directory src/entropy_floor/EAGLE/
and the floor
is applied in the drift and kick operations of the hydro scheme. It is
also used in some of the other subgrid schemes.
The first limit, labelled as Cool
, is typically used to prevent
lowdensity highmetallicity particles to cool below the warm phase because
of overcooling induced by the absence of metal diffusion. This limit plays
only a small role in practice. The second limit, labelled as Jeans
, is
used to prevent the fragmentation of highdensity gas into clumps that
cannot be resolved by the coupled hydro+gravity solver. The two limits are
sketched on the following figure.
An additional overdensity criterion above the mean baryonic density is applied to prevent gas not collapsed into structures from being affected. To be precise, this criterion demands that the floor is applied only if \(\rho_{\rm com} > \Delta_{\rm floor}\bar{\rho_b} = \Delta_{\rm floor} \Omega_b \rho_{\rm crit,0}\), with \(\Delta_{\rm floor}\) specified by the user, \(\rho_{\rm crit,0} = 3H_0/8\pi G\) the critical density at redshift zero [1], and \(\rho_{\rm com}\) the gas comoving density. Typical values for \(\Delta_{\rm floor}\) are of order 10.
The model is governed by 4 parameters for each of the two limits. These are
given in the EAGLEEntropyFloor
section of the YAML file. The parameters
are the Hydrogen number density (in \(cm^{3}\)) and temperature (in
\(K\)) of the anchor point of each floor as well as the powerlaw slope
of each floor and the minimal overdensity required to apply the
limit. Note that, even though the anchor points are given in terms of
temperatures, the slopes are expressed using a powerlaw in terms of
entropy and not in terms of temperature. For a slope of \(\gamma\) in
the parameter file, the temperature as a function of density will be
limited to be above a powerlaw with slope \(\gamma  1\) (as shown on
the figure above). To simplify things, all constants are converted
to the internal system of units upon reading the parameter file.
For a normal EAGLE run, that section of the parameter file reads:
EAGLEEntropyFloor:
Jeans_density_threshold_H_p_cm3: 0.1 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in, expressed in Hydrogen atoms per cm^3.
Jeans_over_density_threshold: 10. # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in.
Jeans_temperature_norm_K: 8000 # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold, expressed in Kelvin.
Jeans_gamma_effective: 1.3333333 # Slope of the EAGLE Jeans limiter entropy floor
Cool_density_threshold_H_p_cm3: 1e5 # Physical density above which the EAGLE Cool limiter entropy floor kicks in, expressed in Hydrogen atoms per cm^3.
Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold, expressed in Kelvin.
Cool_gamma_effective: 1. # Slope of the EAGLE Cool limiter entropy floor
SWIFT will convert the temperature normalisations and Hydrogen number density thresholds into internal energies and densities respectively assuming a neutral gas with primoridal abundance pattern. This implies that the floor may not be exactly at the position given in the YAML file if the gas has different properties. This is especially the case for the temperature limit which will often be lower than the imposed floor by a factor \(\frac{\mu_{\rm neutral}}{\mu_{ionised}} \approx \frac{1.22}{0.59} \approx 2\) due to the different ionisation states of the gas.
Recall that we additionally impose an absolute minium temperature at all densities with a value provided in the SPH section of the parameter file. This minimal temperature is typically set to 100 Kelvin.
Note that the model only makes sense if the Cool
threshold is at a lower
density than the Jeans
threshold.
Chemical tracers¶
The gas particles in the EAGLE model carry metal abundance information in the form of metal mass fractions. We follow explicitly 9 of the 11 elements that Wiersma et al. (2009)b traced in their chemical enrichment model. These are: H, He, C, N, O, Ne, Mg, Si and Fe [2]. We additionally follow the total metal mass fraction (i.e. absolute metallicity) Z. This is typically larger than the sum of the 7 metals that are individually traced since this will also contain the contribution of all the elements that are not individually followed. We note that all of definitions are independent of any definition of solar the solar metallicity \(Z_\odot\) or of any solar abundance pattern.
As part of the diagnostics, we additionally trace the elements coming from the different stellar evolution channels. We store for each particle the total mass coming from all the SNIa that enriched that particle and the metal mass fraction from SNIa. This is the fraction of the total gas mass that is in the form of metals originating from SNIa stars. By construction this fraction will be smaller than the total metal mass fraction. The same tracers exist for the SNII and AGB channels. Finally, we also compute the iron gas fraction from SNIa. This it the fraction of the total gas mass that is made of iron originating from SNIa explosions.
We finally also compute the smoothed versions of the individual element mass fractions, of the total metal mass fractions, and of the iron gas fraction from SNIa.
The chemistry module in src/chemistry/EAGLE/
includes all the arrays
that are added to the particles and the functions used to compute the
smoothed elements.
When a star is formed (see the section Star formation: Schaye+2008 modified for EAGLE below), it inherits all the chemical tracers of its parent gas particle.
In the snapshots, we output for each gas and star particle:
Name  Description  Units  Comments 

ElementAbundance 
Fraction of the gas/star mass
in the different elements

[]  Array of length
9 for each particle

SmoothedElementAbundance 
Fraction of the gas/star mass
in the different elements
smoothed over SPH neighbours

[]  Array of length
9 for each particle

Metallicity 
Fraction of the gas/star mass
in all metals

[]  
SmoothedMetallicity 
Fraction of the gas/star mass
in all metals
smoothed over SPH neighbours

[]  
TotalMassFromSNIa 
Total mass of the gas/star
that was produced by enrichment
from SNIa stars

[U_M]  
MetalMassFracFromSNIa 
Fraction of the total gas/star
mass that is in metals produced
by enrichment from SNIa stars

[]  
TotalMassFromAGB 
Total mass of the gas/star
that was produced by enrichment
from AGB stars

[U_M]  
MetalMassFracFromAGB 
Fraction of the total gas/star
mass that is in metals produced
by enrichment from AGB star

[]  
TotalMassFromSNII 
Total mass of the gas/star
that was produced by enrichment
from SNII stars

[U_M]  
MetalMassFracFromSNII 
Fraction of the gas/star mass
that is in metals produced by
enrichment from SNII stars

[]  
IronMassFracFromSNIa 
Fraction of the total gas/star
mass in iron produced produced
by enrichment from SNIa stars

[]  
SmoothedIronMassFracFromSNIa 
Fraction of the total gas/star
mass in iron produced produced
by enrichment from SNIa stars
smoothed over SPH neighbours

[] 
The stars will lose mass over their lifetime (up to ~45%). The fractions will
remain unchanged but if one is interested in computing an absolute metal mass
(say) for a star, the InitialMass
(see the section
Star formation: Schaye+2008 modified for EAGLE below) of the star must be used.
The chemistry model only requires a small number of parameters to be specified in the EAGLEChemistry section of the YAML file. These are the initial values of the metallicity and element mass fractions. These are then applied at the start of a simulation to all the gas particles. All 9 traced elements have to be specified An example section, for primordial abundances (typical for a cosmological run), is:
EAGLEChemistry:
init_abundance_metal: 0. # Mass fraction in *all* metals
init_abundance_Hydrogen: 0.755 # Mass fraction in Hydrogen
init_abundance_Helium: 0.245 # Mass fraction in Helium
init_abundance_Carbon: 0. # Mass fraction in Carbon
init_abundance_Nitrogen: 0. # Mass fraction in Nitrogen
init_abundance_Oxygen: 0. # Mass fraction in Oxygen
init_abundance_Neon: 0. # Mass fraction in Neon
init_abundance_Magnesium: 0. # Mass fraction in Magnesium
init_abundance_Silicon: 0. # Mass fraction in Silicon
init_abundance_Iron: 0. # Mass fraction in Iron
Whilst one would use the following values for solar abundances (typical for an idealised lowredshift run):
EAGLEChemistry:
init_abundance_metal: 0.014 # Mass fraction in *all* metals
init_abundance_Hydrogen: 0.70649785 # Mass fraction in Hydrogen
init_abundance_Helium: 0.28055534 # Mass fraction in Helium
init_abundance_Carbon: 2.0665436e3 # Mass fraction in Carbon
init_abundance_Nitrogen: 8.3562563e4 # Mass fraction in Nitrogen
init_abundance_Oxygen: 5.4926244e3 # Mass fraction in Oxygen
init_abundance_Neon: 1.4144605e3 # Mass fraction in Neon
init_abundance_Magnesium: 5.907064e4 # Mass fraction in Magnesium
init_abundance_Silicon: 6.825874e4 # Mass fraction in Silicon
init_abundance_Iron: 1.1032152e3 # Mass fraction in Iron
Individual element abundances for each particle can also be read
directly from the ICs. By default these are overwritten in the code by
the values read from the YAML file. However, users can set the
parameter init_abundance_metal
to 1
to make SWIFT ignore the
values provided in the parameter file:
EAGLEChemistry:
init_abundance_metal: 1 # Read the particles' metal mass fractions from the ICs.
The ICs must then contain values for these three fields (same as what is written to the snapshots):
Name  Description  Units  Comments 

ElementAbundance 
Fraction of the gas/star mass
in the different elements

[]  Array of length
9 for each particle

Metallicity 
Fraction of the gas/star mass
in all metals

[]  
IronMassFracFromSNIa 
Fraction of the total gas/star
mass in iron produced produced
by enrichment from SNIa stars

[] 
If these fields are absent, then a value of 0
will be used for all
of them, likely leading to issues in the way the code will run.
Gas cooling: Wiersma+2009a¶
The gas cooling is based on the redshiftdependent tables of Wiersma et
al. (2009)a that include
elementbyelement cooling rates for the 11 elements (H, He, C, N, O,
Ne, Mg, Si, S, Ca and Fe) that dominate the total rates. The tables
assume that the gas is in ionization equilibrium with the cosmic microwave
background (CMB) as well as with the evolving Xray and UV background from
galaxies and quasars described by the model of Haardt & Madau (2001). Note that this model
ignores local sources of ionization, selfshielding and nonequilibrium
cooling/heating. The tables can be obtained from this link
which is a repackaged version of the original tables. The code reading and interpolating the
table is located in the directory src/cooling/EAGLE/
.
The Wiersma tables containing the cooling rates as a function of redshift, Hydrogen number density, Helium fraction (\(X_{He} / (X_{He} + X_{H})\)) and element abundance relative to the solar abundance pattern assumed by the tables (see equation 4 in the original paper). As the particles do not carry the mass fraction of S and Ca, we compute the contribution to the cooling rate of these elements from the abundance of Si. More specifically, we assume that their abundance by mass relative to the table’s solar abundance pattern is the same as the relative abundance of Si (i.e. \([Ca/Si] = 0\) and \([S/Si] = 0\)). Users can optionally modify the ratios used for S and Ca. Note that we use the smoothed abundances of elements for all calculations.
Above the redshift of Hydrogen reionization we use the extra table containing net cooling rates for gas exposed to the CMB and a UV + Xray background at redshift nine truncated above 1 Rydberg. At the redshift or reionization, we additionally inject a fixed userdefined amount of energy per unit mass to all the gas particles.
In addition to the tables we inject extra energy from Helium II reionization using a Gaussian model with a userdefined redshift for the centre, width and total amount of energy injected per unit mass. Additional energy is also injected instantaneously for Hydrogen reionisation to all particles (active and inactive) to make sure the whole Universe reaches the expected temperature quickly (i.e not just via the interaction with the now much stronger UV background).
For noncosmological run, we use the \(z = 0\) table and the interpolation along the redshift dimension then becomes a trivial operation.
The cooling itself is performed using an implicit scheme (see the theory documents) which for small values of the cooling rates is solved explicitly. For larger values we use a bisection scheme. Users can alternatively use a NewtonRaphson method that in some cases runs faster than the bisection method. If the NewtonRaphson method does not converge after a few steps, the code reverts to a bisection scheme, that is guaranteed to converge. The cooling rate is added to the calculated change in energy over time from the other dynamical equations. This is different from other commonly used codes in the literature where the cooling is done instantaneously.
We note that the EAGLE cooling model does not impose any restriction on the particles’ individual timesteps. The cooling takes place over the time span given by the other conditions (e.g the Courant condition).
Finelly, the cooling module also provides a function to compute the temperature of a given gas particle based on its density, internal energy, abundances and the current redshift. This temperature is the one used to compute the cooling rate from the tables and similarly to the cooling rates, they assume that the gas is in collisional equilibrium with the background radiation. The temperatures are, in particular, computed every time a snapshot is written and they are listed for every gas particle:
Name  Description  Units  Comments 

Temperature 
Temperature of the gas as
computed from the tables.

[U_T]  The calculation is performed
using quantities at the last
timestep the particle was active

Note that if one is running without cooling switched on at runtime, the
temperatures can be computed by passing the temperature
runtime flag (see
Command line options). Note that the tables then have to be available as in
the case with cooling switched on.
The cooling model is driven by a small number of parameter files in the EAGLECooling section of the YAML file. These are the reionization parameters, the path to the tables and optionally the modified abundances of Ca and S as well as the flag to attempt using the NewtonRaphson scheme to solve the implicit problem. A valid section of the YAML file looks like:
EAGLECooling:
dir_name: /path/to/the/Wiersma/tables/directory # Absolute or relative path
H_reion_z: 11.5 # Redhift of Hydrogen reionization
H_reion_ev_p_H: 2.0 # Energy injected in eV per Hydrogen atom for Hydrogen reionization.
He_reion_z_centre: 3.5 # Centre of the Gaussian used for Helium reionization
He_reion_z_sigma: 0.5 # Width of the Gaussian used for Helium reionization
He_reion_ev_p_H: 2.0 # Energy injected in eV per Hydrogen atom for Helium II reionization.
And the optional parameters are:
EAGLECooling:
Ca_over_Si_in_solar: 1.0 # (Optional) Value of the Calcium mass abundance ratio to solar in units of the Silicon ratio to solar. Default value: 1.
S_over_Si_in_solar: 1.0 # (Optional) Value of the Sulphur mass abundance ratio to solar in units of the Silicon ratio to solar. Default value: 1.
newton_integration: 0 # (Optional) Set to 1 to use the NewtonRaphson scheme for the explicit cooling problem.
Particle tracers¶
Over the course of the simulation, the gas particles record some information
about their evolution. These are updated for a given particle every time it is
active. The EAGLE tracers module is located in the directory
src/tracers/EAGLE/
.
In the EAGLE model, we trace the maximal tempearature a particle has reached and the time at which this happened. When a star is formed (see the section Star formation: Schaye+2008 modified for EAGLE below), it inherits all the tracer values of its parent gas particle. There are no parameters to the model but two values are added to the snapshots for each gas and star particle:
Name  Description  Units  Comments 

Maximal Temperature 
Mximal temperature reached by
this particle.

[U_T]


Maximal Temperature scalefactor OR
Maximal Temperature time 
Scalefactor (cosmological runs)
or time (noncosmological runs) at
which the maximum value was reached.

[]
OR
[U_t]

Star formation: Schaye+2008 modified for EAGLE¶
The star formation is based on the pressure implementation of Schaye & Dalla Vecchia (2008) with a metaldependent starformation density threshold following the relation derived by Schaye (2004). Above a density threshold \(n^*_{\rm H}\), expressed in number of Hydrogen atoms per (physical) cubic centimeters, the star formation rate is expressed as a pressurelaw \(\dot{m}_* = m_g \times A \times \left( 1 {\rm M_\odot}~{\rm pc^2} \right)^{n} \times \left(\frac{\gamma}{G_{\rm N}}f_gP\right)^{(n1)/2}\), where \(n\) is the exponent of the KennicuttSchmidt relation (typically \(n=1.4\)) and \(A\) is the normalisation of the law (typically \(A=1.515\times10^{4} {\rm M_\odot}~{\rm yr^{1}}~{\rm kpc^{2}}\)). \(m_g\) is the gas particle mass, \(\gamma\) is the adiabatic index, \(f_g\) the gas fraction of the disk and \(P\) the total pressure of the gas including any subgrid turbulent terms.
Once a gas particle has computed its star formation rate, we compute the probability that this particle turns into a star using \(Prob= \min\left(\frac{\dot{m}_*\Delta t}{m_g},1\right)\). We then draw a random number and convert the gas particle into a star or not depending on our luck.
The density threshold itself has a metallicity dependence. We use the smoothed metallicty (metal mass fraction) of the gas (See Chemical tracers) and apply the relation \(n^*_{\rm H} = n_{\rm H,norm}\left(\frac{Z_{\rm smooth}}{Z_0}\right)^{n_{\rm Z}}\), alongside a maximal value. The model is designed such that star formation threshold decreases with increasing metallicity. This relationship with the YAML parameters defining it is shown on the figure below.
In the EAGLE model, the pressure entering the star formation includes pressure from the unresolved turbulence. This is modeled in the form of a polytropic equation of state for the gas \(P = P_{\rm norm}\left(\frac{\rho}{\rho_0}\right)^{\gamma_{\rm eff}}\). For practical reasons, this relation is expressed in term of densities. Note that unlike the entropy floor, this is applied at all densities and not only above a certain threshold. This equation of state with the relevant YAML parameters defining it is shown on the figure below.
To prevent star formation in noncollapsed objects (for instance at high redshift when the whole Universe has a density above the threshold), we apply an overdensity criterion. Only gas with a density larger than a multiple of the critical density for closure can form stars.
Additionally to the pressurelaw corresponding to the KennicuttSchmidt relation described, above, we implement a second density threshold above which the slope of the relationship varies (typically steepens). This is governed by two additional parameters: the density at which the relation changes and the second slope. Finally, we optionally use a maximal density above which any gas particle automatically gets a probability to form a star of 100%.
The code applying this star formation law is located in the directory
src/star_formation/EAGLE/
. To simplify things, all constants are converted
to the internal system of units upon reading the parameter file.
For a normal EAGLE run, that section of the parameter file reads:
# EAGLE star formation parameters
EAGLEStarFormation:
EOS_density_norm_H_p_cm3: 0.1 # Physical density used for the normalisation of the EOS assumed for the starforming gas in Hydrogen atoms per cm^3.
EOS_temperature_norm_K: 8000 # Temperature om the polytropic EOS assumed for starforming gas at the density normalisation in Kelvin.
EOS_gamma_effective: 1.3333333 # Slope the of the polytropic EOS assumed for the starforming gas.
KS_normalisation: 1.515e4 # The normalization of the KennicuttSchmidt law in Msun / kpc^2 / yr.
KS_exponent: 1.4 # The exponent of the KennicuttSchmidt law.
min_over_density: 57.7 # The overdensity above which starformation is allowed.
KS_high_density_threshold_H_p_cm3: 1e3 # Hydrogen number density above which the KennicutSchmidt law changes slope in Hydrogen atoms per cm^3.
KS_high_density_exponent: 2.0 # Slope of the KennicutSchmidt law above the highdensity threshold.
EOS_temperature_margin_dex: 0.5 # (Optional) Logarithm base 10 of the maximal temperature difference above the EOS allowed to form stars.
KS_max_density_threshold_H_p_cm3: 1e5 # (Optional) Hydrogen number density above which a particle gets automatically turned into a star in Hydrogen atoms per cm^3.
threshold_norm_H_p_cm3: 0.1 # Normalisation of the metaldependant density threshold for star formation in Hydrogen atoms per cm^3.
threshold_Z0: 0.002 # Reference metallicity (metal mass fraction) for the metaldependant threshold for star formation.
threshold_slope: 0.64 # Slope of the metaldependant star formation threshold
threshold_max_density_H_p_cm3: 10.0 # Maximal density of the metaldependant density threshold for star formation in Hydrogen atoms per cm^3.
gas_fraction: 1.0 # (Optional) The gas fraction used internally by the model.
Stellar enrichment: Wiersma+2009b¶
Supernova feedback: Dalla Vecchia+2012 & Schaye+2015¶
# EAGLE stellar enrichment and feedback model
EAGLEFeedback:
use_SNII_feedback: 1 # Global switch for SNII thermal (stochastic) feedback.
use_SNIa_feedback: 1 # Global switch for SNIa thermal (continuous) feedback.
use_AGB_enrichment: 1 # Global switch for enrichement from AGB stars.
use_SNII_enrichment: 1 # Global switch for enrichement from SNII stars.
use_SNIa_enrichment: 1 # Global switch for enrichement from SNIa stars.
filename: ./yieldtables/ # Path to the directory containing the EAGLE yield tables.
IMF_min_mass_Msun: 0.1 # Minimal stellar mass considered for the Chabrier IMF in solar masses.
IMF_max_mass_Msun: 100.0 # Maximal stellar mass considered for the Chabrier IMF in solar masses.
SNII_min_mass_Msun: 6.0 # Minimal mass considered for SNII feedback (not SNII enrichment!) in solar masses.
SNII_max_mass_Msun: 100.0 # Maximal mass considered for SNII feedback (not SNII enrichment!) in solar masses.
SNII_wind_delay_Gyr: 0.03 # Time in Gyr between a star's birth and the SNII thermal feedback event.
SNII_delta_T_K: 3.16228e7 # Change in temperature to apply to the gas particle in a SNII thermal feedback event in Kelvin.
SNII_energy_erg: 1.0e51 # Energy of one SNII explosion in ergs.
SNII_energy_fraction_min: 3.0 # Maximal fraction of energy applied in a SNII feedback event.
SNII_energy_fraction_max: 0.3 # Minimal fraction of energy applied in a SNII feedback event.
SNII_energy_fraction_Z_0: 0.0012663729 # Pivot point for the metallicity dependance of the SNII energy fraction (metal mass fraction).
SNII_energy_fraction_n_0_H_p_cm3: 0.67 # Pivot point for the birth density dependance of the SNII energy fraction in cm^3.
SNII_energy_fraction_n_Z: 0.8686 # Powerlaw for the metallicity dependance of the SNII energy fraction.
SNII_energy_fraction_n_n: 0.8686 # Powerlaw for the birth density dependance of the SNII energy fraction.
SNIa_max_mass_Msun: 8.0 # Maximal mass considered for SNIa feedback and enrichment in solar masses.
SNIa_timescale_Gyr: 2.0 # Timescale of the exponential decay of the SNIa rates in Gyr.
SNIa_efficiency_p_Msun: 0.002 # Normalisation of the SNIa rates in inverse solar masses.
SNIa_energy_erg: 1.0e51 # Energy of one SNIa explosion in ergs.
AGB_ejecta_velocity_km_p_s: 10.0 # Velocity of the AGB ejectas in km/s.
SNII_yield_factor_Hydrogen: 1.0 # (Optional) Correction factor to apply to the Hydrogen yield from the SNII channel.
SNII_yield_factor_Helium: 1.0 # (Optional) Correction factor to apply to the Helium yield from the SNII channel.
SNII_yield_factor_Carbon: 0.5 # (Optional) Correction factor to apply to the Carbon yield from the SNII channel.
SNII_yield_factor_Nitrogen: 1.0 # (Optional) Correction factor to apply to the Nitrogen yield from the SNII channel.
SNII_yield_factor_Oxygen: 1.0 # (Optional) Correction factor to apply to the Oxygen yield from the SNII channel.
SNII_yield_factor_Neon: 1.0 # (Optional) Correction factor to apply to the Neon yield from the SNII channel.
SNII_yield_factor_Magnesium: 2.0 # (Optional) Correction factor to apply to the Magnesium yield from the SNII channel.
SNII_yield_factor_Silicon: 1.0 # (Optional) Correction factor to apply to the Silicon yield from the SNII channel.
SNII_yield_factor_Iron: 0.5 # (Optional) Correction factor to apply to the Iron yield from the SNII channel.
Blackhole creation¶
Blackhole accretion¶
Blackhole repositioning¶
AGN feedback¶
[1]  Recall that in a noncosmological run the critical density is set to 0, effectively removing the overdensity constraint of the floors. 
[2]  Wiersma et al. (2009)b originally also followed explicitly Ca and and S. They are omitted in the EAGLE model but, when needed, their abundance with respect to solar is assumed to be the same as the abundance of Si with respect to solar (See the section Gas cooling: Wiersma+2009a) 