Testing geological proxies for deep-time tidal model simulations

Tides are a key driver of a range of Earth system processes, and we now have the capacity to simulate tidal dynamics on a range of temporal and spatial scales. Deep-time tidal model simulations have been used to provide insight into past ocean circulation patterns, evolution of life and the developments of the Earth-Moon system's orbital configuration. However, these tidal model simulations are relatively poorly constrained and validated because of a lack of readily available proxies. The feasibility of using two types of proxy is explored here; (1) sedimentary deposits which can directly estimate palaeotidal ranges, and (2) black shale, to constrain three palaeotidal model simulations for different time slices. Specifically, three palaeotidal range proxies were used for the early Devonian (400 Ma), three palaeotidal range proxies and five black shales for the lower Jurassic (185 Ma), and


| INTRODUCTION
Ocean tides impact a range of Earth system processes.They control the locations of productive shelf sea fronts (Simpson & Hunter, 1974), sustain the climate-regulating global overturning circulation (Wilmes et al., 2021;Wunsch & Ferrari, 2004), drive ocean primary production (Sharples et al., 2007;Tuerena et al., 2019) and set the environment for key evolutionary events (Balbus, 2014;Byrne et al., 2020).The dissipation of tidal energy also slows down Earth's spin and forces the moon to recede to conserve angular momentum (Bills & Ray, 1999;Daher et al., 2021), meaning the tides are a first-order controller of daylength.Recent tidal model results (Green et al., 2017(Green et al., , 2018) ) show significantly less energetic tides in the past.This has far-reaching consequences for the Earth system, for example, its biogeochemical cycles, and may have been a driving force in the oxygenation of the atmosphere (Klatt et al., 2021).
Numerical modelling of palaeotides relies on tectonic reconstructions for boundary conditions (see Green et al., 2022, and references above).However, despite numerous publications outlining characteristics of palaeotides from the palaeobiological and geological records, numerical simulations are poorly constrained as the proxy information is not readily accessible to the modelling community.This is partially rectified by collating information on deeptime tides from different sources and using these data to constrain numerical tidal model simulations for three time slices: 400, 185 and 95 Ma.These were chosen for their representative nature and the availability of suitable proxy data in the literature.The results can also benefit the regional palaeogeographical reconstructions: if the tidal conditions are verified at a location, the regional topography is accurate as well.
Palaeotidal range estimates were obtained from results presented in the literature; palaeotidal range is determined through palaeoenvironmental interpretation, and may be estimated by analogy to modern tidal environments (Klein, 1971;Wells et al., 2005).For instance, small-scale sedimentary structures are usually distributed in mesotidal environments (Reineck, 1975), while large-scale structures or plane beds are in macrotidal settings (Dalrymple et al., 1990).However, these interpretations may only provide a rough estimate of the tidal range, rather than an accurate definition (Collins et al., 2021).Furthermore, well-preserved tidalites display ebb-flood, and spring-neap or spring-neap-like cycles, alongside signals of the diurnal inequality, and can be used to constrain days per month and days per year counts in the geological past (Archer, 1996;Coughenour et al., 2009).Under ideal conditions, tidal range in a specific location can be directly estimated from tidalites where complete, fining-upward sequences of sediment are preserved (Devries Klein, 1971;Slingerland, 1986;Tanavsuu-Milkeviciene & Plink-Bjorklund, 2009;Williams, 2000).For example, the sharp contrast of thin neap couplets and thick spring couplets in tidal bundles suggested macrotidal range or higher (Tessier & Gigot, 1989).In other cases, the palaeotidal range could be estimated from the stratigraphic thickness of intertidal deposits in tidal flat units (Klein, 1970(Klein, , 1971)).Sediment can sometimes also be translated into a tidal current range based on bedload transport rates for non-cohesive sediments (Ward et al., 2020), where fine-grained sediments (silt and clay) are more common at low tidal ranges.However, exploring this aspect falls outside the scope of this study.
Furthermore, tidal mixing fronts separate a vertically mixed water column from a stratified one (Bowers & Simpson, 1987;Simpson & Bowers, 1981;Simpson & Hunter, 1974).Their positions are determined by the tidal current magnitude and the water depth, and they are found near contours of where H is water depth (m), and u is the tidal current magnitude (m 2 /s); χ the commonly referred to as the Simpson-Hunter parameter (m 2 /s 3 ).Consequently, mapping the fronts provides a proxy for palaeotidal current magnitudes.Black shale deposition occurs in poorly ventilated anoxic conditions incompatible with strong tidal currents (Abdi et al., 2021;He et al., 2022;Stow et al., 2001;Wignall & Newton, 2001).It is expected therefore that the presence of a black shale in the geological record will indicate a strongly stratified water column, which can be tracked through the locations of the tidal mixing fronts (Simpson & Hunter, 1974).Note that the specific driver leading to the anoxic event is not important, rather the fact that a well-mixed water column is also well-ventilated and if the water column is mixed, all properties are mixed too.The tides provide a continuous supply of mechanical energy for mixing and a black shale must therefore sit in a stratified regime.A similar method, tracking microfossil assemblages, has also been suggested but will not be pursued here (Scourse et al., 2002).
To date, a few proof-of-concept studies used tidal deposits to constrain numerical tidal model simulations (Byrne et al., 2020;Green et al., 2020;Wells et al., 2010;Zuchuat et al., 2022), but the information used was limited to a few data points for specific regions and time slices.Consequently, this paper presents a systematic comparison between environmental proxies and numerical tidal model results. (1) The simulations of past tides were made using OTISthe Oregon State University Tidal Inversion Software-a dedicated tidal model which has been used extensively to simulate deep-time and present day (PD) tides (Blackledge et al., 2020;Byrne et al., 2020;Egbert et al., 2004;Green et al., 2022).OTIS was benchmarked against other forward tidal models and was shown to perform well (Stammer et al., 2014) which provides a numerical solution to the linearised shallow-water equations forced by the tide only, as presented by the following Equations ( 2) and (3): Here, U = u H is the tidal volume transport (u is the horizontal velocity vector [m/s] and H is the water depth [m]), f the Coriolis parameter (rad/s), g the acceleration due to gravity (m/s 2 ), η the sea-surface elevation (m), η SAL the self-attraction and loading elevation (m), η EQ is the elevation of the equilibrium tide (m) and F the tidal energy dissipation term (W/kg or m 2 /s 3 ).The dissipation is parameterised through two components, denoted F B and F W respectively, representing bed friction and energy losses due to tidal conversion, that is, the energy transferred into a baroclinic tide respectively.Friction is parameterised using the standard quadratic law, F B = C D u|u|, where C D = 0.003 is a dimensionless drag coefficient (Taylor, 1920), whereas the tidal conversion term may be written as F W = C U, with a conversion coefficient, C, expressed by following Equation 4 (Green & Nycander, 2013;Zaron & Egbert, 2006): Here, γ = 50 represents a dimensionless scaling factor representing unresolved bathymetric roughness, N H is the buoyancy frequency at the seabed (unit of s −1 ), N represents the vertical average of the buoyancy frequency (rad/s) and ω the frequency of the tide.The buoyancy frequency, N (rad/s), is given by N 2 = −g/ρ∂ρ/∂z, where ρ is the density (kg/m 3 ).The details of the density field are not known for the period discussed here, so the values of N based on a statistical fit to observed PD values presented by Zaron and Egbert (2006) were used.Consequently, N(x,y) = 0.00524exp(−z/1300), where z is the vertical coordinate (m), and the constants 0.00524 and 1300 have units of s −1 and m respectively.The model results are relatively insensitive to changes in stratification and this parameter space is not explored further (Byrne et al., 2020;Egbert et al., 2004;Green et al., 2020).

| Reconstructions and simulations
The palaeobathymetry data came from Scotese and Wright (2018) and were supplied at 1/10° horizontal resolution in both latitude and longitude.All bathymetries effectively ran from 89° S to 89° N in latitude due to the introduction of land at the poles to handle the convergence of the model grid cells at high latitudes.Note that outside of near-resonant states, tidal simulations are relatively insensitive to small-scale topographic changes and the blocking of the poles (Egbert et al., 2004;Wilmes & Green, 2014).The details for each era are summarised below and described in more detail in each section.
All time slices were simulated for the M 2 , S 2 and K 1 constituents.All time slices were simulated using PD tidal forcing as well as changed forcing to parameters relevant for each time slice (Daher et al., 2021)-see Table 1 for details.The model outputs the amplitudes and phases of the surface elevations, η, and transports, U, for each of the constituents.

| PROXIES
Two types of tidal proxies are explored here: (1) sedimentary deposits that can be used to directly estimate T A B L E 1 Forcing parameters used in the model simulations, with data from Daher et al. (2021).The M 2 forcing factor is based on the change in lunar distance associated with the change in orbital periods.The K 1 forcing factor is made up of 2/3 from the Moon and 1/3 from the Sun.

| Direct proxies: Palaeotidal range
Here, palaeotidal ranges interpreted from tidal deposits in the literature are regarded as 'direct proxies' (DPs); these have been used previously to constrain palaeotidal modelling (Byrne et al., 2020;Green et al., 2020;Zuchuat et al., 2022).Tidal ranges from the literature were categorised into the four standard categories: microtidal (0-2 m), mesotidal (2-4 m), macrotidal (4-6 m) and hypertidal (>6 m; Archer, 2013; note that the oceanographic literature use megatidal for a range >6 m).In this investigation, if the results from the tidal model for a proxy location fall within the category provided by the proxy, the simulation is considered accurate.Palaeotidal range proxies will be used to validate the simulations for 400 Ma and 185 Ma.Details of the direct tidal proxies are summarised below and in Tables S1A and S1B in the Supplementary Information and the locations are presented in Figure 1A,B.Three palaeotidal range proxies for 400 Ma were found in the literature, all also used by Byrne et al. (2020); note that the simulations here use a different reconstruction and higher resolution than Byrne et al. (2020) did.A meso-macrotidal regime was discovered in the fine to coarse-grained, tide-dominated deltaic deposits of the Rezekne and Pärnu formations in the Devonian Baltic Basin; these make up (Direct proxy 1) DP1 (Figure 1A; Tanavsuu-Milkeviciene & Plink-Bjorklund, 2009;Tovmasjana, 2013).Another meso-macrotidal environment was found in the rippled and cross-bedded siliciclastic and dolomitic tidal flat facies of the Padeha Formation in the Tabas Block of the Central-East Iranian Microcontinent (DP2; Wendt et al., 2004;Zand-Moghadam et al., 2014).Griffing et al. (2000) and Rust et al. (1989) inferred a mesotidal regime from the tidally deposited sandstone and mudstone bodies in the Cap-aux-Os Member of the Battery Point Formation in the Catskill Clatic Wedge, Canada (DP3 in Figure 1A).For the 185 Ma time slice, three direct proxies were found in the literature (Figure 1B).Sellwood (1972) determined the minimum tidal range as 1 m from the thickness of sandstone channel-fill sequences in Gry's Lower Coal Series in Bornhom, Denmark (DP4).A macrotidal regime was concluded based on the dimensions of estuarine bedforms in the incised valleys of the Ostreaelv Formation (DP5) in the Niell Klinter Group in Jameson Land, Greenland (Ahokas et al., 2014).Lastly, the Helsingborg Member of the Gassum Formation in southern Sweden and the Galgeløkke Member of the Rønne Formation on Bornholm (both combined as DP6) consisting of tidal flat and channel facies, were deposited in a microtidal to mesotidal environment (Nielsen et al., 1989).

| Indirect proxies: Black shale
It has been suggested that tidal rhythmites, which can indicate palaeotidal range, are predominantly formed within middle to inner estuaries (Tessier, 2023).This can introduce uncertainties when validating tidal model results because the reconstruction is unlikely to cover small-scale estuaries.Furthermore, where direct proxies have been found for the tidal regimes, these are almost all meso-macrotidal, leaving low tidal ranges without proxies (at this stage).To address this, it is proposed that locations of black shale constitute an additional tidal proxy for tides in a shelf sea setting.
Black shale is an indirect proxy which can ultimately constrain tidal current velocities.Tide-driven mixing controls stratification shelf seas, with a tidal mixing front marking the point between vertically mixed and stratified areas as discussed in the introduction.Mapping front positions from the model output and comparing them to locations with black shale therefore constitutes a validation metric: the shale must sit on the stratified side of the front and if they do not, the model simulation is probably incorrect.The identified palaeo-locations of the black shale formations used here are shown as 'BS' in Figure 1B,C.A total of 13 locations were identified: five for the 185 Ma timeslice and a further eight for 95 Ma.Note that black shale deposits in deep water will not be considered here, as anoxia in the deep ocean is not controlled by tides.Furthermore, black shale deposited between the Precambrian (400 Ga) to the Devonian (419-358 Ma) was deposited at a time when ocean chemistry was significantly different and there is evidence to suggest that there was not enough oxygen in the marine environment to form oxic waters (Aharon, 2005;Kimura & Watanabe, 2001).

| PD model validation
The performance of the PD set-up of the tidal model was evaluated compared to the TPXO9 satellite altimetry constrained product (Egbert & Erofeeva, 2002).The predicted M 2 and S 2 tidal amplitudes for PD are presented in Figure 2, with the globally spatial-averaged root-meansquare (RMS) errors for M 2 and S 2 amplitudes calculated to be 9.8 cm and 4.4 cm respectively.

| 400 Ma
The predicted M 2 and S 2 amplitudes and the mean spring tidal range at 400 Ma are presented in Figure 3A,B.A region with high amplitudes of M 2 and S 2 is situated in the western and northern parts of Laurussia, southern Siberia and the north-east of Gondwana, where the M 2 amplitude exceeds 1.5 m, while S 2 amplitude exceeds 1 m.Less energetic waters are found in east Laurussia and northern Siberia as M 2 or S 2 microtidal regimes dominate these areas.The model prediction here is generally consistent with the lower-resolution simulations for 400 Ma from Byrne et al. (2020), but variations in tidal predictions occur in specific areas due to differences in the utilised palaeobathymetry data.For instance, this study reveals

F I G U R E 6
The M 2 tidal range indicated by direct proxies and the corresponding model prediction for 185 Ma.Note that in this figure, DP4 has been moved from the original location on land to the nearest coastal ocean grid cell.significantly higher 2 amplitudes in southern Siberia (up to 2.5 m) and lower S 2 amplitudes in north-east Laurussia (lower than 0.5 m) compared to the prediction from Byrne et al. (2020).
The mean spring tide range is computed as 2 M 2 + S 2 , as shown in Figure 3C,D, where M 2 and S 2 are the corresponding tidal amplitudes for principal lunar semidiurnal (M 2 ) and principal solar semidiurnal (S 2 ).Macrotidal areas (4-6 m) are located along the west and north coastline of Laurussia, southern Siberia and north-east region of Gondwana.In contrast, mesotidal regions (2-4 m range) were found on the south-east coast of Laurussia.Compared with direct tidal proxies collected for 400 Ma and plotted in Figure 4, the model prediction matches reasonably well with the proxies at DP3 and DP2 (see Table S1A for details).However, the simulation does not agree with the tidal proxy of the DP1.DP1 was located in a meso-macrotidal delta, whereas the simulation indicates a microtidal environment in that region.This discrepancy may be an effect of model resolution: the Pärnu and Rēzekne formations were deposited in transitional fluvial-tide-dominated flats, tidal channels or in a deltaic distributary channel that, at the current model resolution, is not resolved.

| 185 Ma
The 185 Ma model results are validated using both tidal range proxies and black shale (see Table S1B for a summary of proxy information).It should be noted that DP4 (Gry's Lower Coal Series) will not work as a tidal proxy because it is too far inland.This is of course still useful information because it tells us that the reconstruction needs to be modified to encompass the proxy location.These simulations host a relatively weak global tide, except for a macrotidal regime in the western Tethys Sea (Figure 5).It can be argued that if DP4 is moved to the nearest coastal grid cell, it fits well with the tidal model result of 1.4 m M 2 tidal range (Figure 6).The predicted tidal range DP6 (Helsingborg Member the Gassum Formation and the Galgeløkke Member of the Rønne Formation) is 1.25 m, which is in the range of the tidal proxy there (Figure 6).However, the tidal range prediction at the DP5 site (Ostreaelv Formation) is 0.30 m, which is a considerable underestimate compared to the 4-6 m macrotidal range was identified in the proxy.It is proposed that this is another resolution issue with the model grid: a macrotidal range in a shallow seaway may not be fully resolved with a 1/10° model resolution, resulting in the omission of a resonant feature in the inner part of the seaway.The calculated Simpson-Hunter criterion for 185 Ma in Figure 7 and Table 2 shows that the model result matches four of five black shale tidal proxies (BS2-BS5).BS1 is located at a borderline mixed region with a value of log 10 χ ~ 2 in the model simulations, whereas the proxy of course points to a stratified water column.This could be remedied by increasing the depth of the area to reduce current speeds and increase χ, as demonstrated by Equation (1).Again, the tidal proxies give information about the quality of the tectonic reconstructions as well as acting as a validation tool for the tides.

| 95 Ma
The predicted M 2 amplitude and its corresponding velocity magnitude for the 95 Ma time slice are presented in Figure 8.In general, this is a quiescent time slice, with weak tides in the vast epeiric seas covering PD Africa, Asia and Europe.The exceptions are high-velocity zones northeast of PD Madagascar and north Australia (Figure 8A).
There is a lack of direct tidal proxies for 95 Ma and instead the Simpson-Hunter criterion and black shale formations will be used as proxies.The results shown in Figure 9 indicate that much of the 95 Ma shelf seas were stratified and that the black shale records are matched by the modelled stratification in six of eight locations (see Table 3).The two mismatched locations, BS10 and BS12, are again located at the boundary of mixed and stratified regions with values of log 10 χ = 2.0 and 1.8 respectively.A small correction of the depth would also ensure a stratified water column at these two locations.The positive correlation between black shale palaeo-locations and tidal front locations suggests that black shale can serve as indirect proxies for palaeotides.

| DISCUSSION
Simulations of the tides for three deep-time time slices are presented with the results validated with two types of geological tidal proxies: palaeotidal ranges deducted from tidal deposits (a direct proxy) and black shale (an indirect proxy).Direct proxies were collected from published literature for 400 and 185 Ma; for both time slices, the model performed reasonably well with proxies agreeing in 2/3 of the locations.However, in areas where it fell short, it is argued that the bathymetry of the reconstruction could be modified to ensure a better fit.Using direct proxies is the preferred method for validating palaeotidal models (Byrne et al., 2020;Zuchuat et al., 2022) but they are rare in the literature.Therefore, it is proposed to use black shale as an indirect proxy, providing an upper limit of the tidal current speeds, with a proof-of-concept study presented for time slices at 185 and 95 Ma.Because of the interconnections between tidal current speeds, stratification and potential for anoxia, it is argued that in cases where the Simpson-Hunter parameter denotes a stratified water column (i.e.log 10  > 2-2.3), the presence of black shale in that region can be attributed to the water-mass stratification.The results and proxies agree in 10 of 13 locations across both time slices.It would be easy to change F I G U R E 8 (A) Simulated M 2 tidal amplitudes; (B) simulated M 2 tidal current magnitude for the 95 Ma time slice.For clarity, the proxy locations are marked in Figure 9.
the water depth the model and agree.This way, both verified tidal model simulations and improved reconstructions are obtained.
The main uncertainty in this type of work is in the bathymetric reconstructions because the tidal dynamics are largely controlled by the bathymetry (Byrne et al., 2020;Green et al., 2017Green et al., , 2020;;Zuchuat et al., 2022).The uncertainty or error of the model simulations is given by the RMS value provided; this also highlights that the main uncertainty is the bathymetry.It is very difficult to quantify the uncertainty or accuracy of the proxies, because they are proxies and the modern analogues are largely missing.However, this is largely mitigated by the range of the tidal characteristics from the proxies and if the model simulation falls within that range, there can be confidence in the results for both the model and the proxy.The problem is when the two do not agree and at least one of the two-the model or the proxy-is incorrect.There is no way of knowing which one at this stage, and more work is needed to improve the model set-up, for example, by using higher resolution in the model simulations and constraining the reconstructions better.This work is underway and left for a future publication.

| CONCLUSION
While it can be argued that the results make sense from a dynamical perspective, the encouraging correlations between the model results and proxies show that the model is reasonably correct, and that the methodology works.It also demonstrates that there could be a wealth of viable tidal proxies available in the literature, and that collecting and collating them is a worthy effort to constrain deep-time tidal simulations.kinds of proxies explored here, and it is worth investigating further potential proxies found in the literature.For example, grain size distributions could be used alongside bedforms, such as ripples, to provide direct constraints on the current speed (Baas, 1999;Baas et al., 2021;Davis & Dalrymple, 2012;Oost & Baas, 1994).Other proxies can come from palaeobiology, where species distributions can tell us about the size of intertidal zones (and hence tidal range) and, again, help constrain the bathymetry.Matching information of geological formations and basins can also contain estimates of shelf width and specific topography which can be used to further reconstruct bathymetries.The same is true for palaeocurrent directions which could potentially be used to verify shoreline trends.These investigations are left for future publications.

ACKNO WLE DGE MENTS
Simulations were done on Supercomputing Wales with funding from HEFCW.BG and JAMG acknowledge funding from the Natural Environment Research Council (MATCH, NE/S009566/1).LAMF, JMH and OP were recipients of a 2021 Bangor University Undergraduate internship; they contributed equally to the manuscript and are listed in alphabetical order.We express our sincere gratitude for the constructive and insightful comments provided by Valentin Zuchuat, an anonymous reviewer and the editors.

F
I G U R E 1 The palaeogeographical reconstructions and the location of tidal proxies for (A) 400 Ma; (B) 185 Ma; (C) 95 Ma, with the specifics of direct proxies (DP) and black shales (BS).F I G U R E 2 (A) The simulated M 2 tidal amplitudes in metres for the PD control simulations; (B) The simulated S 2 tidal amplitudes.

F
I G U R E 3 (A) Simulated M 2 tidal amplitudes; (B) simulated S2 tidal amplitudes; (C) simulated mean spring tidal range for 400 Ma, calculated by 2 M 2 + S 2 , and the marked palaeotidal range proxies; (D) close-up of ocean region surrounding Laurussia.F I G U R E 4 The tidal range indicated by direct proxies and the corresponding model prediction for 400 Ma.The modelled tidal range is the range in the gridcell nearest to the proxy location, where the error bar shows the largest and smallest values in a 3 × 3 grid box centred on the proxy location.

F
I G U R E 5 (A) Simulated mean spring tidal amplitudes for the 185 Ma time slice; (B) close-up of the Laurasian Seaway where the proxies are located.

F
The model prediction of Simpson-Hunter criterion () and associated stratification state for the 185 Ma time slice.

F
Predicted Simpson-Hunter criterion for 95 Ma, and the palaeo-location of black shale proxies.T A B L E 3 The model prediction of Simpson-Hunter criterion (logarithms to 10) and the associated tidal stratification for 95 Ma.The locations of black shales are expected to sit in a stratified water column.