Seismic records contain information that allows geoscientists to make inferences about the structure and properties of the Earth’s interior. Traditionally, seismic imaging and tomography methods require wavefields to be generated and recorded by identifiable sources and receivers, and use these directly-recorded signals to create models of the Earth’s subsurface. However, in recent years the method of seismic interferometry has revolutionised earthquake seismology by allowing unrecorded signals between pairs of receivers, pairs of sources, and source-receiver pairs to be constructed as Green’s functions using either cross-correlation, convolution or deconvolution of wavefields. In all of these formulations, seismic energy is recorded and emitted by surrounding boundaries of receivers and sources, which need not be active and impulsive but may even constitute continuous, naturally-occurring seismic ambient noise. In the first part of this thesis, I provide a comprehensive overview of seismic interferometry, its background theory, and examples of its application. I then test the theory and evaluate the effects of approximations that are commonly made when the interferometric formulae are applied to real datasets. Since errors resulting from some approximations can be subtle, these tests must be performed using almost error-free synthetic data produced with an exact waveform modelling method. To make such tests challenging the method and associated code must be applicable to multiply-scattering media. I developed such a modelling code specifically for interferometric tests and applications. Since virtually no errors are introduced into the results from modelling, any difference between the true and interferometric waveforms can safely be attributed to specific origins in interferometric theory. I show that this is not possible when using other, previously available methods: for example, the errors introduced into waveforms synthesised by finite-difference methods due to the modelling method itself, are larger than the errors incurred due to some (still significant) interferometric approximations; hence that modelling method can not be used to test these commonly-applied approximations. I then discuss the ability of interferometry to redatum seismic energy in both space and time, allowing virtual seismograms to be constructed at new locations where receivers may not have been present at the time of occurrence of the associated seismic source. I present the first successful application of this method to real datasets at multiple length scales. Although the results are restricted to limited bandwidths, this study demonstrates that the technique is a powerful tool in seismologists’ arsenal, paving the way for a new type of ‘retrospective’ seismology where sensors may be installed at any desired location at any time, and recordings of seismic events occurring at any other time can be constructed retrospectively – even long after their energy has dissipated. Within crustal seismology, a very common application of seismic interferometry is ambient-noise tomography (ANT). ANT is an Earth imaging method which makes use of inter-station Green’s functions constructed from cross-correlation of seismic ambient noise records. It is particularly useful in seismically quiescent areas where traditional tomography methods that rely on local earthquake sources would fail to produce interpretable results due to the lack of available data. Once constructed, interferometric Green’s functions can be analysed using standard waveform analysis techniques, and inverted for subsurface structure using more or less traditional imaging methods. In the second part of this thesis, I discuss the development and implementation of a fully non-linear inversion method which I use to perform Love-wave ANT across the British Isles. Full non-linearity is achieved by allowing both raypaths and model parametrisation to vary freely during inversion in Bayesian, Markov chain Monte Carlo tomography, the first time that this has been attempted. Since the inversion produces not only one, but a large ensemble of models, all of which fit the data to within the noise level, statistical moments of different order such as the mean or average model, or the standard deviation of seismic velocity structures across the ensemble, may be calculated: while the ensemble average map provides a smooth representation of the velocity field, a measure of model uncertainty can be obtained from the standard deviation map. In a number of real-data and synthetic examples, I show that the combination of variable raypaths and model parametrisation is key to the emergence of previously-unobserved, loop-like uncertainty topologies in the standard deviation maps. These uncertainty loops surround low- or high-velocity anomalies. They indicate that, while the velocity of each anomaly may be fairly well reconstructed, its exact location and size tend to remain uncertain; loops parametrise this location uncertainty, and hence constitute a fully non-linearised, Bayesian measure of spatial resolution. The uncertainty in anomaly location is shown to be due mainly to the location of the raypaths that were used to constrain the anomaly also only being known approximately. The emergence of loops is therefore related to the variation in raypaths with velocity structure, and hence to 2nd and higher order wave-physics. Thus, loops can only be observed using non-linear inversion methods such as the one described herein, explaining why these topologies have never been observed previously. I then present the results of fully non-linearised Love-wave group-velocity tomography of the British Isles in different frequency bands. At all of the analysed periods, the group-velocity maps show a good correlation with known geology of the region, and also robustly detect novel features. The shear-velocity structure with depth across the Irish Sea sedimentary basin is then investigated by inverting the Love-wave group-velocity maps, again fully non-linearly using Markov chain Monte Carlo inversion, showing an approximate depth to basement of 5 km. Finally, I discuss the advantages and current limitations of the fully non-linear tomography method implemented in this project, and provide guidelines and suggestions for its improvement.