A Nineteen Day Earth Tide Measurement with a MEMS Gravimeter

The measurement of tiny variations in local gravity enables the observation of subterranean features. Gravimeters have historically been extremely expensive instruments, but usable gravity measurements have recently been conducted using MEMS (microelectromechanical systems) sensors. Such sensors are cheap to produce, since they rely on the same fabrication techniques used to produce mobile phone accelerometers. A significant challenge in the development of MEMS gravimeters is maintaining stability over long time periods, which is essential for long term monitoring applications. A standard way to demonstrate gravimeter stability and sensitivity is to measure the periodic elastic distortion of the Earth due to tidal forces - the Earth tides. Here we present a nineteen day measurement of the Earth tides, with a correlation coefficient to the theoretical signal of 0.979. The estimated bias instability of the proposed gravimeter is 8.18 microGal for an averaging time of ~400 s when considering the raw, uncompensated data. The bias instability extracted from the sensor electronic noise sits just under 2 mircoGal for an averaging time of ~200 s. After removing the long-term temperature and control electronics effects from the raw data, a linear drift of 268 microGal/day is observed in the data, which is among one of the best reported for a MEMS device. These results demonstrate that this MEMS gravimeter is capable of conducting long-therm time-lapse gravimetry, a functionality essential for applications such as volcanology.


Introduction
Gravimetry has been used extensively over the past century in several fields. Due to the high cost of commercial instruments (around £70k for a portable device) the oil and gas industry has been the most prolific user of gravimetry, where it is often used as a pre-drilling survey method to investigate subterranean geological features 1;2 . Gravimetry has also been used, however, for sinkhole analysis 3 , finding tunnels and cavities for the defence sector 4;5 , CO 2 sequestration 6 , geothermal reservoir monitoring 7 , archaeology 8 , hydrology 9;10 , and volcanology 11;12;13;14;15 .
The last item in this list is perhaps the field in which gravimetry could offer the most societal benefit, but for which its use has been severely limited by the high capital cost of the equipment. Gravimetry is the only means by which mass variations within volcanoes can be measured.
Gravimetry therefore offers a great advantage to hazard forecasting because mass changes often precede eruptive processes 16;12;17 . Freire et al. 18 state that 'more than 8% of the world's 2015 population lived within 100 km of a volcano with at least one significant eruption'. The significance to humanity of improving eruption forecasting cannot therefore be overstated. Gravimetry has been used to a limited extent in this setting, but spatial resolution has been limited by the number of instruments that can be affordably bought. Whilst single sensors can only provide a point measurement, if the price of gravimeters could be significantly reduced then multi-pixel gravity imaging would become possible. An analogy could be drawn between this and the example of a digital camera; with a single pixel you cannot capture a picture, you only have a light sensor, but this soon changes as you gain pixels. In 2019 Carbone et al. 19  Due to the inverse square law of gravity, a single measurement does not have a single unique solution; it is impossible to tell whether a signal has a source that is big and far away, or small and close. Inversion techniques therefore need to be applied to understand the data. Arrays such as the planned network on Mt. Etna present a great opportunity to reduce the problem of inversion, because by using two or more sensors one benefits from triangulation in reducing the number of possible solutions. The more sensors one introduces, the more information that it is possible to glean from the data. In essence, the data from an array of sensors is greater than the sum of its parts; and multi-pixel gravity imaging has the potential to change the way that mass changes within volcanoes are monitored.

Design of the MEMS Gravimeter
The MEMS gravimeter has two distinct components -a monolithically etched silicon sensor with a high acceleration sensitivity; and an integrated capacitive readout scheme to measure the sensor's output. Here, the design aspects of the core MEMS sensor will be discsused briefly, followed by an explanation of the implemented readout scheme.
The design of the sensor is a based upon the device first presented by the group in 2016 22;23;24 .
Significant changes have been made since this time, however, in order to miniaturise the sensor and improve both its sensitivity and stability. The sensor is a relative gravimeter: it is able to measure changes in gravity by monitoring the relative displacement of a proof-mass on a spring. For an oscillating system such as this, the input acceleration, a, and the corresponding proof-mass displacement, z, are directly proportional to each other through the relationship a = 4π 2 f 2 z, where f is the fundamental resonant frequency, and the bandwidth of the sensor. Here, the mass and the supporting springs are etched monolithicically from a single piece of silicon using standard photolithography and etching techniques 25 (see figure 1). The four symmetric springs utilise a geometric anti-spring design 26;27;28;29;30 . The details of how this design can be tuned for the purposes of gravimeter fabrication is detailed in the paper by Middlemiss et. al. 31 . This design enables a decreased spring stiffness in the vertical direction whilst limiting the motion in the horizontal and out-of-plane axes. Decreasing the stiffness -and hence the oscillation frequency -in the operation axis is important because it means that a proof-mass will move more for a given change in gravity. This is particularly useful in cases where the sensitivity of the system is limited by the displacement readout. Softening the spring, however, comes at a cost because the robustness of the sensor decreases as the resonant frequency drops. With the intention of developing a much more robust sensor, the spring design was optimised to obtain a working resonant frequency of 7.35 Hz (compared to 2.2 Hz in earlier publications 22 ). Given that there was also a desire to improve on the sensitivity reported in the earlier work, it was necessary to improve the displacement sensitivity measurement of the proof-mass.
The measurement of the MEMS displacement is, of course, another fundamental function of the device. In Middlemiss et. al. 22 , an optical shadow sensor was used to measure the displacement of the mass 32 . An LED was used to cast light on a photodiode, and the MEMS was placed in the light beam. As the mass moved, the shadow cast on the photodiode altered the photocurrent, which was used as a proxy for displacement. This methodology was limiting for several reasons. A displacement sensitivity of less that 1 nm was difficult to achieve, ultimately limiting the acceleration sensitivity of the overall sensor. The components needed to be mounted on a bespoke machined block of fused silica to reduce temperature sensitivity. The disadvantage of this block was twofold: it was expensive to machine, and it was far too big to fit inside a standard MEMS package. Whilst this device configuration was used to conduct field tests 33;34;35 , further minaturisation was required.
To facilitate the device miniaturisation, a capacitive displacement method was adopted to readout the motion of the proof-mass. A set of metal comb electrodes were patterned on the top surface of the proof-mass while a complementary set of electrodes were patterned on a fixed glass layer that was separated from the proof-mass layer by a gap of 40 µm using SU8 spacers. The combs on the two layers formed overlapping capacitors. While the proof-mass comb was electrically driven by a 40 kHz differential pair of sine waves, the signal was picked up by the set of fixed combs on the glass layer. Any displacement of the proof-mass modulates the overlapping capacitance and, thus, the output current flowing through the capacitors. Therefore, by monitoring the changing overlapping capacitance, the position of the mass could be measured. The current obtained at the output of the overlapping capacitors was fed to an operational amplifier that was configured in the conventional transimpedance amplification (TIA) topology. The transimpedance topology, despite its simplicity, is very effective in reading current from high impedance sources because its signal gain is insensitive to capacitive parasitics at the pick-up stage * . The voltage signal obtained after the amplification stage was then fed to a bench-top lock-in amplifier. The signal was demodulated * The noise gain of the TIA stage, however, is a function of the input parasitics and needs more attention during circuit design when the sensitivity of the system is limited by the electronics noise floor. and low-pass filtered, and the signal amplitude that carried the displacement information was acquired and logged onto a PC through a Labview routine at a sampling frequency of roughly 0.2 Hz. Using the modulation-demodulation lock-in approach to extract the signal amplitude helps in lowering the noise floor as it only allows the noise at the carrier frequency to pass through while rejecting the noise at all other frequencies.
The capacitive displacement measurement mechanism is advantageous both because it offers a significant reduction in size compared to the shadow sensor (both the sensor and the readout were implemented on the same chip stack), but also because it had a displacement sensitivity that was significantly better than the optical shadow sensor. With a capacity to resolve 50-100 pm of proof-mass displacement, the capacitive displacement readout represented a sensitivity improvement over the shadow sensor by a factor of 10-20. This readout methodology could be improved even further because the displacement sensitivity of the capacitive readout ( δC δz , where C is the overlapping capacitance) is a function of the comb geometry as well as that of the gap between the two overlapping combs 36 . By increasing the comb finger density and bringing the proof-mass and the glass pick-up layer closer, the displacement sensitivity can be enhanced. In the case of our sensor, however, we opted for design parameters (comb density, and gap) that provided a significant safety margin on the fabrication and assembly tolerances. This, of course, can be optimised in the future.
To ensure the long-term stability and the noise performance of the MEMS gravimeter, it was necessary to control the temperature of the device. As discussed in Middlemiss et. al. 22 , tempera- 3 Results

Earth Tide Measurement
One means of assessing the performance of a gravimeter is to measure the Earth tides. The

Spectral Analysis
The spectral analysis of the gravimeter data provides further insights into the presence of periodic signals of varying frequencies in the acceleration data. The blue series of figure 3 is a presentation of the drift-corrected, regressed data (the light grey series of figure 2) in the form of an amplitude spectral density graph. To complement the slow-sampled gravimeter data, a second series of Figure 2 : Two time series plots of the Earth tide measurement, conducted at Glasgow University in April 2019. In both graphs, the red series is a theoretical signal calculated using T-Soft, the gray-scale series are the experimental data recorded using the MEMS gravimeter. From lightest grey to black, these series represent increasing filtering of the data using a low-pass filter (LPF). The black series was filtered with a 20 µHz bandwidth LPF. All of the theoretical data is presented after a regression analysis was used to remove drifts caused by temperature variations within the system. The upper graph shows the full data set, the lower graph shows a zoomed in selection of this data.
higher-frequency gravimeter data (sampled at a rate of 20 Hz) is also included in the plot † . This can be observed in the orange series in figure 2. The yellow series in figure 2 is the electronic noise floor of the sensor.
The blue series in figure 2 has a clearly identifiable double peak at 10 −5 Hz. These two peaks represent the diurnal, and semi-diurnal periodicity of the Earth tides. This signal splitting occurs as the relative phase of the Sun-Earth-Moon system changes. When all are aligned in a line, the gravitational signal of the Sun and the Moon acts in phase on the Earth. When these three bodies form a right-angled triangle, however, the gravitational signal of the Sun and Moon acting on the Earth are out of phase. It is this same phase drift that causes the total amplitude to vary between spring and neap tides. Another demonstration of this Earth tide signal splitting has not before been presented in the literature (to the knowledge of the authors). Such data has not previously been presented because its observation requires the a clean time series of well over a week. In addition to the double-peak, there is a small spurious peak occurring around 25 mHz. This signal is not geophysical in nature, and has been identified as being caused by the temperature control electronics. The roll-off seen in the slow-sampled data at higher frequencies is due to the low-pass filter operation performed by the lock-in amplifier.
In the orange series, further signals can be observed. The furthest right of these is the resonance peak of the MEMS mass-on-spring system at 7.35 Hz. This represents the maximum frequency at which the device can record meaningful data. The sensitivity to signals of higher frequencies than this resonance will rapidly decrease due to the nature of a simple harmonic oscillator. The second set of signals that can be identified in the orange series is geophysical in nature: the primary and secondary mircoseisms. The primary microseism occurs at frequencies between 0.04 and 0.15 Hz 41 , † The gravimeter data was not sampled at a higher sampling rate for the entire bandwidth to avoid unnecessary accumulation of data. and is caused by interactions between ocean waves and the sea bed 42 . The secondary microseism has a larger amplitude, and occurs between 0.08 and 0.3 Hz 41 . This signal is also related to ocean activity, but is caused by interactions in shallow water between wind-driven ocean waves, and waves reflected back from the coastline. Both the primary microseisms vary in amplitude and phase, depending on ocean conditions.
As the higher-frequency band of the sensor output is dominated by the ground motion, it is difficult to estimate the true noise floor of the sensor. One way to get rid of this issue is to use a reference seismometer to effectively subtract the sesimic signals in the high-frequency band.
However, we did not have access to such an instrument. We instead estimated the sensitivity by measuring the true electronic noise floor of the full sensor set-up. This was achieved by grounding the modulation sine wave drive at the input of the MEMS capacitors. The measured sensor noise data is represented by the yellow series in the figure, and shows a flat spectral response for the measured frequency bandwidth. The roll-off observed for higher frequencies, again is due to the low-pass filtering. The sensitivity of the gravimeter is estimated to be 18 µGal/ √ Hz at 1 Hz, but at this frequency the sensor is limited by seismic noise. As discussed in the next section, the ultimate noise floor of the sensor is 0.91 µGal, which can be achieved over an integration time of 250 s.
It should also be noted that the device is also a sensitive seismometer. One way in which his can be observed is via the presence of anthropogenic noise in 2. It is quite evident that the amplitude of the high-frequency components increases and then decreases with the changing footfall within the building. Since it is difficult to capture this time-localised information using the standard fourier analysis, we have used wavelet analysis to isolate the periodic variation in the anthropoegenic noise. The details of the analysis are included in the Appendix. Figure 3 : An amplitude spectral density plot of the post-regression instrument data (the blue and the orange series), and the instrument noise data (the yellow series). The blue series represents data that were recorded at a sampling rate of 0.18 Hz. The orange and the yellow series were recorded at a sampling rate of 20 Hz. Three distinct frequency bands of interest can be observed in the data -these are highlighted using the grey bands. At 10 −5 Hz, the diurnal and semi-diurnal components of the Earth tides are visible; the primary and secondary microseisms can be observed at around 0.2 Hz; and the resonant peak of the MEMS device itself is located at 7.35 Hz. The sensor noise has a flat spectral response for the measured bandwidth. Note: The downward roll-off at high frequencies for the blue and the yellow series is due to post-demodulation low-pass filtering.

Allan Deviation Analysis
To understand the stability and the noise characteristics of the gravimeter, the Allan-Deviation (AD) values (σ) were computed from the raw, uncompensated gravimeter data. ‡ The σ values for the slow-and the fast-sampled data, and the sensor noise are plotted in Figure 4. We have used the fast-sampling data (the orange series in figure 3) for evaluating the stability of the gravimeter for shorter sampling/averaging periods (τ ), and the slow-sampling data (the blue series in figure   3) for assessing the long-term stability. In the case of the slow-sampled data, the first few τ values (τ = 5.5 s, 11 s) have been excluded because the low-pass filter was shaping the frequency response.
As the fast-sampled data was logged for a sufficiently long duration, there were enough samples in the fast-sampled data to compute the σ values for the missing τ values in the slow-sampled series.  Dashed series are a result of fitting straight lines across different averaging periods for each data set. The slope of each fit is represented in the form of power-laws. Note: The small integration time (τ ) data is not shown for the blue (τ < 2 2 t slow , where, t slow = 5.5 sec.) and the yellow series (τ < 2 3 t fast , where, t fast = 0.05 sec.) as the AD values for these times are affected by the low-pass filtering stage.

Drift characteristics
As it is difficult to visualise the hidden features in the the drift of the device by simply computing the AD values, we used regression analysis and polynomial curve fitting on the raw, uncompensated data to study the drift behaviour. The raw sensor data is plotted in Figure 5(a) and clearly shows a non-linear behaviour (the rounded 'M' shape) on top of a linear trend. To remove the tides from affecting the analysis, we applied low-pass filtering, with a 2 µHz cut-off, to the data. Setting this bandwidth was a compromise between avoiding filtering the non-linear aspects of the drift while attenuating the tide signal enough from the series. The effect of filtering can be seen in Figure 5(b), where the low-pass filtered data is represented by the blue series. In the next step, we regressed the filtered data against the instrument's temperature and PID control output channels.
Correcting for the long term temperature effects resulted in the removal of the non-linear trend leaving an almost linear drift (the red series) in the data. We next fitted a straight line to the corrected data and obtained a linear drift of 268µGal/day with an r 2 coefficient of 0.99. § It is believed that the origin of this trend is environmental (for example, temperature gradient effects that are not captured through the temperature sensors), and not geophysical.

Conclusions
In this paper, a 19 day measurement of the Earth tides has been presented. The correlation coefficient between this measurement and a theoretical signal is 0.96, demonstrating a remarkable stability for a MEMS sensor. This sensor utilises an in-plane capacitive displacement readout, which has reduced the size of the device compared to previous work by the authors. This displacement sensor has also been used to increase the displacement sensitivity to 40-50pm, leading § The drift rate obtained using this analysis is slightly higher compared to the one obtained using the AD approach, and could be attributed to less than ideal signal processing (regression, fitting) of the data.

Figure 5 :
The raw data (a) and the drift analysis (b) plots for the gravimeter. The raw data is first low-pass filtered to retain only very long-term features (<2 µHz) in the raw data (blue series in (b)). The extracted drift is then processed to regress out the impact of the temperature and the control instrumentation. The post-regression data (red series) is then fitted against a straight line (dashed black series) to obtain the drift rate of the gravimeter.
to a bias instability value of 8.18µGal. A network of these sensors is currently being constructed for deployment on Mount Etna as part of the NEWTON-g collaboration 20 , to create the first multi-pixel gravity imager.
• P. R. Utting conducted the cross correlation analysis of the tidal time-series.
• G. D. Hammond had oversight of the experimental design, led the initial conceptual design of the gravimeter, conducted data analysis, and contributed to the manuscript preparation.

Author information
• The research data relevant to this paper will be stored on the University of Glasgow's Enlighten Repository on publication • The authors have no competing financial interests.
• Correspondence and requests for material should be addressed to abhinav.prasad@glasgow.ac.uk, richard.middlemiss@glasgow.ac.uk, or giles.hammond@glasgow.ac.uk