Thermally-driven scintillator flow in the SNO+ neutrino detector

This paper will report numerical simulations of laminar convective flow within a spherical container, carried out using OpenFoam to investigate a flow phenomenon observed within the SNO+ liquid scintillator particle detector. It complements several earlier papers [1, 2, 3] addressing a similar subject, and although the novelty of the work lies in its context rather than its methodology, the simplicity of the problem (in terms of domain geometry, fluid homogeneity and thermal forcing) implies that the results we present should have some generality, supplementing what is a surprisingly sparse literature [4] on a common and important type of flow [5, 6] occurring in laboratory beakers, storage tanks and so forth.

For present purposes there is no need to delve deeply into the SNO+ science goals and methodology (for details see [7, 8, 9]). The SNO+ detector is located in a rock cavity some 2 km underground near Sudbury, Ontario. The scintillator fluid, a linear alkylbenzene (LAB), is contained within an acrylic sphere of 6 m radius (Fig. 1), and rises up a cylindrical ‘neck’ (radius 0.75 m) to the interface with a ‘cover gas’ of ultrapure N₂ within the Universal Interface (UI), about 13 m above the centre of the sphere. The detector sphere (hereafter ‘AV’ for ‘acrylic vessel’) and neck ‘float’ within a surrounding mass of purified water that fills the rock cavity to a level about 1 m below the top of the neck, hold-down ropes compensating for the buoyancy of the AV (the scintillator’s density is nominally about $860\,\mathrm{kg\,m^{-3}}$). Within the water, approximately concentric with the AV and at a radius of about 8.5 m from its centre, some 9400 inward-looking photomultipler tubes (PMTs), mounted on a frame referred to as the PSUP (PMT Support), respond to photons emerging from the detector — and carry the information from which physics data are deduced.

It is characteristic of this type of detector that one seeks to discriminate a relatively small number of target events, e.g. scattering of an electron by a solar neutrino, or a hypothetical neutrinoless double beta decay, from a ‘background’ of vastly more numerous events that stem mostly from natural radioactivity within and near the detector. Accordingly every effort is made to minimize radioactive contamination, by optimally selecting material parts, and by filtration and purification of the scintillator and external water. Nevertheless in practice the scintillator will carry a low level of activity, particularly from short-lived progeny of the ${}^{238}\mathrm{U}$ and ${}^{232}\mathrm{Th}$ chains during and after filling and recirculation operations, and this is monitored by quantifying the associated event rates. The contaminant of primary interest in this paper is radon (${}^{222}\mathrm{Rn}$, half-life 3.8 days), which is present in the laboratory air at a level of about $10^{2}\;\mathrm{Bq\,m^{-3}}$ [7]. Within the AV the concentration of ${}^{222}\mathrm{Rn}$ can be quantified from the decay rates of its daughter radionuclides, i.e. the transformation of the daughter bismuth (${}^{214}\mathrm{Bi}$) by beta decay to polonium (${}^{214}\mathrm{Po}$) and the subsequent alpha decay of the ${}^{214}\mathrm{Po}$, a sequence producing a delayed coincidence signal that is easily identified. It is critical for SNO+ that the radon level in the detector be many orders of magnitude below that in ambient air or water, which is achieved by isolating the scintillator from room air to the maximum possible degree¹¹1At present the ratio of the radon concentration in the cover gas to that in the air is smaller than $2\times 10^{-4}$ [7], and further improvement is anticipated. Please note that the level of BiPo-214 activity evident in Figs. (2,5) below was seen just after the completion of scintillator operations to mix in the fluor (2,5-diphenyhloxazole or ‘PPO’, incorporated at a concentration of 2.2 g/L as a wavelength-shifter), and is much higher than, and not characteristic of, the level now achieved..

Evidently then, in the context of the purity of liquid scintillator particle detectors — SNO+ being one of several in operation — one is going to be concerned to understand and perhaps control scintillator motion, which may transport a contaminant into the “fiducial volume”, i.e. a subregion of scintillator sufficiently distant from the wall as to exclude background events occurring in or near the AV walls. Absent any pumping or intentional agitation, scintillator motion can be driven only by thermal gradients or temporal trends, and in practice such convective motion fields are very weak. In this regard one may take the example of the Borexino (liquid scintillator) underground neutrino detector, which is similar in size and type to SNO+, and whose scintillator motion patterns have been investigated in a sequence of papers [10, 1, 2]. Early measurements gave proof of thermally-driven motion within the detector, driven largely by seasonal variation of laboratory air temperature, and modifications were made to reduce the amplitude of that motion (insulation of the outer wall of the detector, and provision of active temperature controls) with the result that lateral inhomogeity of the detector’s boundary temperature was reduced to a level of order 0.1 K. Even so, the scientific goals of the Borexino experiment necessitate quantification of the weak convective motion²²2These authors note the paucity of existing studies of the “natural convection problem under consideration, represented by fluid flow inside [a] closed, stratified, near-equilibrium system”., as reported by [2].

Returning to the SNO+ experiment, it had been envisaged that other than during intervals of ‘AV recirculation’ (i.e. extraction, purification and reintroduction of LAB), stable thermal stratification of scintillator in the neck would mitigate convective transport from the gas/liquid interface at the top of the neck down into the fiducial volume. In that context, and motivated by observation of an unanticipated pattern of scintillator motion, it is the purpose of this paper to present numerical simulations that illustrate and typify the types of flow that may occur in the detector in response to idealized but plausible thermal inhomogeneities.

The outline of the paper is as follows. Section 2 will detail an actual contamination event, a quiescent period (i.e. no pumping) during which a transient and decaying radon layer was observed to slowly sink within the AV. Section 3 will briefly cover OpenFoam (‘OF’, the open source software used to simulate flow in the detector) and its method of application. Direct scintillator velocity measurements are unavailable, and so in Section 4 simulations of a closely analogous flow, viz. that within a 3 L spherical flask of water subjected to controlled warming from its boundary, is compared with corresponding measurements [11]. Confidence in the methodology having been established, Section 5 will show that modest wall heat fluxes (of order $0.1\,\mathrm{W\,m^{-2}}$) induce buoyancy-driven motion of scintillator within the spherical AV that explains (semi-quantitatively) what has been observed, and that if such forcing combines with stable temperature stratification of the scintillator fluid a high frequency wave motion ensues. Finally Section 6 focuses on mixing down the neck, as driven by distinct types of thermal disturbance, and examines the impact of bulk thermal stratification upon that mixing.

In early June 2021, just after the LAB fill had been completed, an interesting event was observed in the pattern of BiPo-214 decays within the SNO+ detector. After some 5 h of AV recirculation on 31 May, a radon-contaminated layer had formed³³3The ${}^{222}\mathrm{Rn}$ entered with the AV recirculation operations, as expected and as previously seen by other scintillator experiments. During ‘normal’ recirculation, scintillator is extracted from the bottom of the AV and re-injected at the base of the neck. Experience has shown that the behaviour of the re-introduced scintillator depends on its temperature relative to that within the AV. Newly returned fluid may promptly undergo mixing, or may form a transient blob or layer at the top of the AV sphere that later sinks., rather uniform in its activity (as seen in the pattern formed by daily-summation of BiPo-214 decays) and spanning approximately the layer $2.5\leq z\leq 5$ m (the coordinate origin used here is the centre of the AV sphere, and $z$ is the vertical coordinate). During the subsequent two weeks, until the next interval of AV recirculation on 14 June, the AV was subjected to no disturbance other than the nitrogen ($\mathrm{N}_{2}$) bubblers⁴⁴4One bubbler emits at the base of the AV, and another at the base of the neck. The bubble stream consists of small – circa. 1 cm diameter – single bubbles, with a cumulative volume flow rate (at atmospheric pressure) of about 150 L in two weeks., and any unmeasured thermal inhomogeneity or trend.

It was observed that over those two weeks the contaminated layer slowly sank, decaying as expected owing to the short (3.8 day) half-life of ${}^{222}\mathrm{Rn}$, but more or less retaining its layered form. Fig. (2) plots that slow descent, as captured by daily integrations of the Bi-214 and Po-214 decays plotted in $(z,\rho^{2})$ coordinates, where $\rho=\sqrt{x^{2}+y^{2}}$.

The radon ‘blob’ of 31 May - 14 June 2021 was not the first such to be observed, but it stands out because of its occurrence during a long interval without external disturbance (the normal situation was for recirculation operations to be paused only on weekends). The need to understand the origin of these events and their evolution prompted the present work, a numerical study of thermally-driven circulation in the AV. What is sought is an initial thermal state, and a choice of thermal forcing on the AV wall, that results in the (assumed) static initial state developing a convective circulation that achieves the observed ‘translation’ of the blob.

Calculations were performed using OpenFoam (‘OF’) v2006⁵⁵5From OpenCFD (www.openfoam.com), an affiliate of Engineering Systems International (ESI)., running under Linux (OpenSuSe Leap 15.2 in single processor mode, or Ubuntu 20.04 in multi-processor mode). The surface shell for each chosen flow domain was generated using Salome⁶⁶6From Open Cascade, www.salome-platform.org to produce ‘stl’ (surface triangle language) files, from which OF’s blockMesh and snappyHexMesh utilities generated the three-dimensional computational mesh⁷⁷7The mesh cells are polyhedra, the majority (more specifically) hexahedra.. Far from domain boundaries the cell size was uniform, and controlled by a setting in OF’s ‘blockMeshDict’ file. Near the walls, layers of cells were added (following a prescription set in ‘snappyHexMeshDict’) to enhance resolution.

The simulations described here used OF’s solver for time-dependent compressible flows, ‘buoyantPimpleFoam’, with an added field variable $C$ representing ${}^{222}\mathrm{Rn}$ and obeying the conservation equation

where $\mathbf{U}$ is the velocity vector; $\rho$ is the fluid density, modelled⁸⁸8In this approximation density is independent of pressure. Some simulations, not reported here, used the ‘incompressible’ solver ‘buoyantBoussinesqPimpleFoam’ (of OF v2006), albeit with an advection term added to the energy equation to parameterize the influence of an implicit uniform and unvarying ‘background’ temperature stratification. as

where $\beta_{T}$ is the thermal expansion coefficient; $\mathcal{D}_{c}$ is the molecular diffusivity of $C$ (prescribed as $\mathcal{D}_{c}=\nu/\mathrm{Sc}$, Sc being the Schmidt number); and $\tau$ is the decay lifetime, in the case of ${}^{222}\mathrm{Rn}$ equal to $4.76\times 10^{5}$ s. Tested in an unforced simulation whose strongly stable stratification mitigated against evolution away from a static initial state, the Rn layer, without moving, simply decayed with the imposed half-life.

Apart from the added field variable $C$, the only other modification made to the solver was to incorporate the Coriolis effect. All simulations assumed laminar flow, which amounts to an assumption that the flows were adequately resolved, both spatially (on the given mesh) and temporally (with the time step applied), obviating any necessity to account for subgrid motion. Further details will be provided below, where they are specific to a given simulation and/or not necessarily obvious even to readers who may be familiar with OpenFoam.

Chow and Akins [11] reported an experiment in which water, contained in a 3 L spherical flask⁹⁹9The neck of the flask was inclined at 45^∘ to the vertical. For the calculations here, the neck was neglected and the radius $R$ of the sphere computed (from the 3 L volume) to be $R=8.95$ cm., was subjected to controlled heating. The flask was immersed within a continuously-stirred water bath, the entire system initially being isothermal at $5^{\circ}$ C. The temperature of the water bath was then progressively increased in such a manner as to sustain a time-invariant temperature difference $T_{\mathrm{b}}-T_{\mathrm{c}}$ between the water bath and the centre of the flask. Suspended in the flask water were glass spheres, whose displacements were measured by a photographic system arranged to provide information on a vertical plane through the centre of the sphere. The authors reported that this fluid system evolves into a pseudo-steady state (i.e. the motion is steady, though the temperature field steadily warms), and presented a single transect of the vertical velocity of the water along an equatorial radius corresponding to the condition that $\Delta T=\overline{T}_{\mathrm{w}}-T_{\mathrm{c}}=2.5$ K, where $\overline{T}_{\mathrm{w}}$ denotes the mean temperature of the inner wall, computed by the authors from measurements and known physical properties¹⁰¹⁰10It is possible that the $2.5^{\circ}\mathrm{C}$ temperature difference had been intended by the authors to refer to the control variable $T_{\mathrm{b}}-T_{\mathrm{c}}$, however certainly they (Chow & Akins) define $\Delta T$ as “fluid temperature difference between inside wall and center of sphere”, and they state in the caption for their Fig. 8 that $\Delta T=2.5^{\circ}\mathrm{C}$.. The value of the Rayleigh number

corresponding to the measured velocity transect was $\mathrm{Ra}=6.5\times 10^{6}$ ($D_{T}$ being the thermal diffusivity; other variables as defined above), and from variations of their experimental configuration Chow and Akins reported that the flow in the flask was laminar provided $\mathrm{Ra}<10^{7}$.

Chow & Akins interpreted (and perhaps designed) their experiment in the context of a prior numerical study [12] of free convection in a sphere, and their measurement provides here the opportunity to gain confidence in the application of OF, in the context of a flow similar in its geometry and forcing to the case of the SNO+ detector: simulations of the laboratory beaker experiment and of the SNO+ neutrino detector differ only as regards the length scale of the mesh, and the physical properties of the fluid. Regarding the latter, Table  (1) specifies the OF thermophysical model used in all the simulations of this paper. Temperature-dependent properties (density, specific heat, viscosity, Prandtl no.) are represented by polynomials, and for the laboratory beaker simulations, the medium being water, polynomials provided by [13] (HolzmannCFD) were used.

The domain boundary for simulations of the Chow-Akins flow consisted of a single spherical ‘patch’ (arbitrarily named ‘shell’), and the desired initial and boundary conditions for temperature were imposed in the OF file ‘caseFolder/0/T’ per the specification of Table (2). By way of explanation, folder ‘0’ at top level within the ‘case folder’ contains a file specific to each dependent variable, and that file prescribes the initial and boundary conditions in standard OF parlance. It is worth noting that, because the simulation holds the entire surface of the domain boundary at $T_{\mathrm{c}}+2.5^{\circ}\mathrm{C}$, it cannot exactly parallel the experiment, for in the latter the flow itself (or rather, the experimental setup in its entirety) determined the inner wall temperature, which could not have been exactly isothermal owing to the anisotropy imposed by gravity.

Fig. (3) compares the measured Chow-Akins velocity transect with two simulations:

• 1.

  Simulation ‘lores’ was performed on a mesh having a total of 309,290 cells, including those within three layers of cells impressed (by the OF utility snappyHexMesh) to represent the near wall region. A 1st-order accurate discretization (‘Gauss upwind’) was used for convective terms, i.e. terms of form $\nabla\cdot(\mathbf{U}\phi)$ where $\mathbf{U}$ is the velocity, with components (UX,UY,UZ) along axes (X,Y,Z) in OF notation, and $\phi$ stands for the convected property.

• 2.

  Simulation ‘hires’ used a mesh with 779,506 cells, there being 10 cell layers impressed to cover the wall region, the outermost of which had a depth about 0.01 mm (i.e. $\sim R/10^{4}$). For this simulation a 2nd-order accurate discretization was used for convective terms.

Overall Fig. (3) shows that the agreement of the two solutions with each other and with the measurements is very good. On their graph, Chow & Akins extrapolated from their outermost (largest $R$) measurement back towards zero, implying they had not measured velocities closer to the wall than that outermost datum. Thus it need not be thought that there is a discrepancy between the simulations and the data between the outermost measured velocity and the wall. As is evident from Fig. (4), showing the simulated velocity profile near the domain wall, certainly the higher-resolution OF run provided good resolution of the near wall flow, and it compares nicely with a solution (for this same problem) given by [5], who numerically integrated the non-dimensional vorticity and temperature equations in polar coordinates, pre-supposing azimuthal symmetry. From Fig. (4) one sees that the boundary-layer Reynolds number $\mathrm{Re}=Ud/\nu$ (where $d\sim 0.5$ mm is the wall boundary layer depth, $U\sim 2\;\mathrm{mm\,s^{-1}}$ the ‘free stream’ velocity and $\nu\sim 10^{-6}\;\mathrm{m^{2}\,s^{-1}}$ the kinematic viscosity) is $\mathrm{O}[1]$, compatible with laminar boundary-layer motion.

Given the near equivalence of the OF model and mesh configurations needed to simulate this small scale experiment and the SNO+ detector, the results of this section provide reassurance that OF should be a satisfactory tool to examine thermally-driven SNO+ flows¹¹¹¹11Two provisos must be listed here: (i) that an adequate computational mesh be provided, and (ii) that the distinction in length scale (between SNO+ and the Chow-Akins experiment), viz. a factor of $6/0.09=67$, is not so large as to imply that the regimes of flow in the two cases (under the given type of forcing) are qualitatively dissimilar. According to the simulations, boundary layer Reynolds numbers are of the order of unity for both systems.. In that regard, and anticipating later results, it is interesting to note that (according to the simulation of the Chow-Akins experiment) the warm wall boundary layer does not ‘detach’ from the wall over the lower hemisphere, and similarly in computations of a ‘reverse Chow-Akins flow’ (i.e. with the wall being held colder rather than warmer than the centre of the flask) the cool wall boundary layer was found to remain attached to the wall over the upper hemisphere.

Simulations of this section relate to the sinking radon layer of Sec. (2), with the domain of the calculations encompassing the entire AV, including the neck. The rate of sink (of the 24-h averaged pattern) was of the order of 0.25 m/day, and the mechanism seems likely to have been the transfer of fluid volume within the wall boundary layer from beneath to above the contamination layer. It is of interest, then, to determine an initial state and forcing that will replicate what had been observed. As no direct measurements of motion in the SNO+ AV exist, and the true initial state and thermal boundary conditions are inaccessible, the purpose of the following simulations is not to gain a highly accurate solution for a specific (but unmeasured) case, but rather, to outline what qualitative patterns of motion exist – and their qualitative effect as demonstrated by the displacement and mixing of a contamination layer from its initial position.

The puzzle presented by the sinking radon layer (Fig. 2) was the absence of any apparent transport mechanism, i.e. in particular the level base of the sinking layer seemed to suggest either that the 24h-averaged pattern was hiding an eddy motion that was accomplishing the transport, or, that fluid volume was being invisibly transferred across/through the contaminant layer along the AV axis or walls. A clue as to the nature of the flow forcing was provided by the observation that in early June 2021, i.e. during the phenomenon, the LAB/cover-gas interface was rising. This could be explained only by thermal expansion of the LAB and demanded an effective heat addition rate of about 300 W, a figure that implies a mean heat flux density over the AV wall of the order of $1\;\mathrm{W\,m^{-2}}$.

Table (3) lists the values adopted for the material properties of LAB. Assuming the motion patterns of interest will primarily be thermally-driven, the property that stands out as being of most importance is the thermal expansion coefficient ($\beta_{T}$). The most uncertain of the properties listed is the Schmidt number (ratio of kinematic viscosity to mass diffusivity), used to evaluate the mass diffusivity $\mathcal{D}_{\mathrm{c}}$ of a species $C$ whose transport equation (Eq. 3) was added to the solver to represent movement and decay of ${}^{222}\mathrm{Rn}$. This however has no bearing on the computed pattern of flow, and unless the Peclet number $U\ell/\mathcal{D}_{\mathrm{c}}$ were spectacularly small, convective transport by the bulk velocity field will be overwhelmingly more imortant than diffusion. The value adopted for the kinematic viscosity $\nu$ is nominal (note: $U$ and $\ell$ are characteristic velocity and length scales).

Within the folder ‘casefolder/constant’ a file (typically named ‘thermophysicalProperties’, but for some solvers named ‘transportProperties’) provides the numerical values of these material properties, along with the user’s choices for the equation of state and the thermodynamic energy variable to be adopted (enthalpy or sensible energy). Table (4) gives the ‘mixture’ block of the OF ‘thermophysicalProperties’ dictionary, for the SNO+ simulations to be shown.

Unless stated otherwise the simulations to follow were made on what will be termed the ‘medium mesh,’ comprising a total of 661,536 cells, and including three layers of cells adjacent to the walls. The depth of the outermost layer of cells was 2.7 cm, and gridlengths based on minimum, mean and maximum cell volume ($V$) were $h=V^{1/3}=(0.009,0.11,0.28)\;\mathrm{m}$. The radon contamination layer initially spanned $2.5\leq z\leq 5$ m, with $C=1$ within the layer and $C=0$ elsewhere. The wall boundary condition for $C$, ‘zeroGradient’, assured there would be no flux to or from walls. (Note: $C$ is named ‘Conc’ on some of the figures to follow.)

The AV neck being the only path, aside from scintillator operations, by which radon is observed to enter the AV sphere, it is of interest to focus on flow and transport over that limited volume of the detector, i.e. a cylinder of radius 0.75 m and spanning $6\leq z\leq 12.75$ m. Here it is relevant to clarify the thermal conditions prevailing outside the neck, as governed by its immersion in water over approximately the lower 6 m of its length but in cover gas over an uppermost section of about 1 m in length. The cover gas temperature is approximately that of the laboratory air, about $17-20^{\circ}$C, whereas the water is controlled at a considerably lower temperature of about $12^{\circ}$C by “cavity recirculation”. The later process extracts water from the top of the column and returns it, after chilling, through two streams that respectively re-enter within and outside the PSUP structure, cooling water bodies referred to as “PSUP water” and “cavity water”. It had been expected that this warm uppermost layer lying atop much cooler fluid would suppress convective motion in the neck, thereby limiting the rate of contaminant ingress and ensuring there would be at least a factor of 50 reduction in the ${}^{222}\mathrm{Rn}$- level in the scintillator at the bottom of the neck relative to the top.

Simulations of this section cover two simple and plausible forms of destabilizing thermal non-uniformity that induce convective mixing down the neck; more complex forms of thermal forcing are easily imagined, but those covered below suffice. The liquid/gas interface at $z=12.75$ m is treated as a ‘free slip/no leak’ surface, and the sidewall and bottom boundary as ‘no slip/no leak’ (i.e. $\mathbf{U}=0$).

Apart from numerous geophysical and planetary studies, the literature on convective flow of a materially homogeneous fluid within a sphere is surprisingly sparse. The present study of motion in the SNO+ liquid scintillator neutrino detector has been motivated by what had been a puzzling phenomenon, and was initiated with the goal of providing a plausible explanation for that motion. The numerical method (OF’s standard solver ‘buoyantPimpleFoam’) being well tried, we have relied on the demonstrated consistency of simulations with measurements in a geometrically congruent small scale experiment (Section 4) as testimony to the realism – qualitative or better – of the full scale (SNO+) calculations.

The simulations suggest that the phenomenon prompting the study, i.e. the slow descent within the detector of a radon contamination layer, was driven by volume exchange from beneath to above that layer, the needed upward volume flux ($Q_{\mathrm{vol}}^{+}$) taking place due to and within a thin convective boundary layer on the AV wall. Consistency of the simulation with the observed multi-day phenomenon hinges to some extent on the absence of any firm constraint in specifying the forcing for the simulated flow, but our choice, viz. a weak heat flux density $q=0.1\,\mathrm{W\,m^{-2}}$ into the fluid from the wall of the upper hemisphere, is not dissimilar from what is implied by (sporadically observed) thermal expansion of the scintillator fluid. Interestingly too, short episodes of contamination rising and spreading laterally from the base of the detector, seen after intervals of so-called ‘reverse’ circulation of the scintillator fluid – ‘reverse’ meaning the LAB is reintroduced at the base of the detector, rather than the base of the neck – are consistent with OF simulations that impose a cooling heat flux around the lower hemisphere of the detector.

The azimuthal symmetry of the SNO+ detector may suggest to readers that the full 3-dimensional (3D) treatment adopted here was uneconomic, and certainly theoretical studies of convection within a sphere (e.g. [5, 12, 14]) have presupposed azimuthal symmetry. Over the course of the work some early simulations adopted as the domain just a single quadrant of the detector, but the additional planar boundaries made mesh generation a more complex task; and as to reducing the problem to 2D, although this would be fine in the context of steady state treatment, in 2D one loses temporal fidelity. In that regard, it appears that there are two or more timescales for the development of this flow. The shortest timescale $\tau_{\mathrm{bl}}$ ($\sim 1\;\mathrm{hr}$, or smaller) relates to the establishment of motion in the wall boundary layer, and so soon as total time $t\gg\tau_{bl}$ (i.e. a few hours after initiation) the sink rate of the blob is established, albeit decreasing with increasing $t$ as the blob sinks, with associated increase of its planar surface area. On a longer time scale ($\tau_{\mathrm{bulk}}$) gross thermal stratification (in response to the wall heat flux) and quasi-horizontal swirl (presumably originating in the Coriois term) are established in the bulk of the AV.

Overall this work indicates that even what might seem to be ‘small’ thermal inhomogeneities or trends can suffice to bring about significant motion (up to $\mathcal{O}[\mathrm{mm/sec}]$) in this type of fluid ‘body’, and the fact of these simulations being for a spherical geometry probably does not limit this suggestion to that particular shape. The motion engendered by some forms of thermal disturbance can be mitigated by ensuring that a stabilizing vertical thermal gradient is sustained, but in this regard the duration and strength of the forcing disturbance are important – a sustained buoyancy flux is liable, eventually, to induce complete mixing. In short, if scintillator motion is to be avoided then (horizontal) thermal homogeneity must be assured to a very fine level.