The thin blue line: A review of shoreline dynamics across time scales and environments

A major thread of theoretical research on the response of shorelines to changing boundary conditions has adapted the moving‐boundary approach from heat transfer and solidification/melting. On sufficiently short time scales, shorelines respond to changes in relative sea level in a simple, geometrically predictable way. On longer time scales, their behaviour becomes far more complex and interesting, because the surface over which the shoreline moves is itself continually modified by morphodynamics that depend sensitively on shoreline location. This makes the shoreline the archetype of moving‐boundary problems in morphodynamics, and subject to potentially counterintuitive behaviours over time scales on which the sediment surface modifies itself as relative sea level changes. We review existing moving‐boundary theories and propose two significant extensions to allow inclusion of first‐order effects of waves and tides. First, we show how wave effects can be included via planform diffusion linked to river‐mouth location, which results in shoreline smoothing during delta‐lobe growth and localized transgression after channel abandonment. Tides produce a low‐gradient region in which the sea and land overlap; we show how this can be treated in a moving‐boundary framework by replacing the shoreline with a ‘mushy region' so that the handoff from land to water occurs over a zone rather than a line. We also propose that the moving‐boundary approach can be readily generalized to other dynamic moving boundaries, such as those separating different regimes of river transport. The shoreline thus serves as a prototype for modelling dynamic facies boundaries along the whole source–sink system. © 2019 The Authors. Earth Surface Processes and Landforms published by John Wiley & Sons Ltd.


Introduction
From a human perspective, the shoreline is perhaps the most striking natural boundary to be found on Earth, separating the familiar terrestrial realm from an alien aqueous realm that in some ways is less well known than the surfaces of other planets. From the standpoint of surface dynamics, the boundary represented by the shoreline is equally striking: it separates a realm of processes we can easily observe directly, like rivers, soil movement and aeolian transport, from one of sea-bed processes associated with waves and currents that are much harder to observe directly ( Figure 1). Yet, in the framework of planetary time, this most dramatic of process boundaries is also one of the most ephemeral. The recent history of glacial advance and retreat implies excursions of shoreline across more or less the whole width of continental shelves over time scales of the order of 10 5 years. The earliest studies of the sedimentary archive, including the celebrated observations of Leonardo Da Vinci, pointed to major changes in sea level and shoreline throughout Earth's history. Yet, as we illustrate in this review, the dynamics of this crucial boundary can be remarkably subtle. Before outlining the main topics for this review, we first explore the origin of complexity in the relation between shoreline and sea level.
Anyone who has spent time on a beach subject to changing tides has a sense of the simple geometric relation between sealevel change and shoreline location: the distance of shoreline change ds is set by the sea-level change dH SL and the beach (terrestrial) slope S t : The simplicity and apparently obvious nature of this relation applies not only to sandcastles, but also to inundation maps that show coastal areas that would be flooded by a given rise in sea level. Given current worries about global warming and its effects on sea level, it is understandable that such maps are proliferating in the mainstream media as well as among specialists (Blum and Roberts, 2009). Our experience at the beach, however, is misleading. Just as we cannot see far enough out to sea to appreciate that the Earth is not flat, we also cannot see far enough back in time to appreciate how the simple geometric view of the shorelinethe 'inclined-plane model'might be wrong. Before explaining, we remind the reader of the distinction between eustatic and relative sea level. Eustatic sea level is the global mean sea level. Relative sea level adds in the effects of local changes in land-surface elevation caused mainly by crustal subsidence or uplift and subsurface compaction. We might then think of rescuing Eq. (1) by replacing H SL with H rsl , where the latter includes these effects. While certainly an improvement, two fundamental problems remain: first, the response to relative sea-level change is also affected by the rate of local erosion or deposition, which certainly depends on the motion of the shoreline inasmuch as the latter divides two transport regimes that could hardly be more different; and second, to the extent that the slope S t is formed of sediment, it too should be influenced by shoreline migration. Since the problem is indeed coupled in these two important ways, it is perhaps better to rethink the premise of Eq. (1) entirely. Doing exactly this is the focus of a line of research that has been going on since 1978, which we will review and build on in this paper. In 1978, Walter Pitman (1978 published a strikingly original analysis of shoreline dynamics which suggested that on planetary time scales, shoreline tracks not eustatic sea level but rather its time rate of change. Shoreline response to cyclic sea-level change is thus 90°out of phase with sea level. Pitman's focus was shoreline response on long time scales. The complete behaviour of his modelincluding the short time limitwas worked out by Angevine (1989). Angevine showed that the inclined-plane model [Eq. (1)] is recovered as the time scale of sea-level change decreases, and that in the long-period limit the out-of-phase behaviour is accompanied by damping of the response amplitude, so that the coastal system acts as a high-pass filter to sea-level cycles. The final result for shoreline migration s(t) over (undifferentiated) coastal-plain and continental shelf topography of slope S, driven by sinusoidally varying sea level H SL = a SL sin (ω SL t) is where θ = ω SL T R and T R is a reference time scale that can be interpreted as the time required for the basin to reach a steady-state balance between subsidence and sedimentation under constant forcing. Pitman's model is essentially geometric, like Eq. (1), but the key conceptual advance is that the sediment surface is not fixed, but rather moves vertically as the crust subsides, sediment is redistributed and sea level changes. The main assumptions are (1) passive margin-style rotational subsidence, (2) constant slope over the coastal plain-continental shelf surface, (3) a sedimentation pattern that maintains this slope against the rotational subsidence and (4) constant sedimentation rate at the shoreline. Pitman originally assumed zero sedimentation and erosion at the shoreline, but Angevine (1989) showed that relaxing this to allow for constant sedimentation at the shoreline did not change the main results of the model. Apart from the phase shift and damped response behaviour, the other key result of Pitman's and Angevine's work was the introduction of a characteristic basin-response time scale for the handoff from 'inclined plane'-type geometric response to the damped, phase-shifted response. This time scale is essentially the same as the equilibrium time defined by Paola et al. (1992) as the time required for the transport surface to achieve a steady state in response to an imposed subsidence profile. We refer the interested reader to the references cited above for further details of the model formulation and its implications. What is important here is not even so much the main resultthe high pass-filtered, phase-shifted shoreline response to eustatic cyclesbut rather that it was, to our knowledge, the first serious attempt to question the apparently obvious geometric link between sea-level change and shoreline migration. The importance of a model like this has little to do with whether its predictions turned out to be correct (we discuss comparisons with field and laboratory studies below). Rather, the real value of Pitman's work is that it changed the way we think about the shoreline. It freed the research community from a conceptual straitjacket, rooted in ordinary experience but incorrect when extended to planetary time and space scales, and opened the way to other models of shoreline dynamics in which the shoreline location is found as part of the solution for the evolution of the sediment surface, and in which the shoreline appears as a dynamic boundary that divides two regimes with quite different transport characteristics. The general approach is termed the method of moving boundaries. In this paper we provide an analytical review of moving-boundary shoreline models, and attempts to test them; propose two new methods that extend their scope; and finally argue that dynamic boundaries over the whole source-sink system can be treated via the same approach with minor modification making all such boundaries, in effect, metaphorical shorelines.

Shoreline Models and Moving Boundaries
Background An early advance in modelling shoreline as a dynamic process on long time scales was a model of prograding deltas developed by Kenyon and Turcotte (1985). The basic idea is that the fluvial transport on the delta topset and the mass flow-dominated transport on the delta foreset can both be modelled as diffusional processes, but with different diffusion coefficients. The shoreline thus provides the location where an inboard higher (fluvial) diffusion coefficient hands off to an outboard lower (mass-flow) coefficient, reflecting the steeper slopes required to transport sediment in the submarine realm, without the help of rivers. The evolution of the boundary emerges from solution of the two diffusion equations, with their respective boundary conditions, along with a coupling condition at the boundary to ensure that mass is conserved and the topography is continuous.
The type of coupling condition used at the shoreline in the Kenyon and Turcotte (1985) model is, though not presented as such, an example of a moving boundary, a concept well known in other areas of analysis such as solidification theory (Crank, 1984). A moving boundary is a boundary within the domain of interest across which there is an abrupt change in the dynamics, as for example at the boundary between a liquid and a solid. Conventional boundary conditions for the equations describing the dynamics within each domain are applied at the moving boundary, along with continuity conditions (e.g. for temperature or elevation). When all conditions are accounted for, the system is formally overconstrained and the extra constraint is then used to solve for the location of the moving boundary at each instant in time.
The moving-boundary method for shorelines The first formulation of shoreline dynamics as a moving-boundary problem is that of Swenson et al. (2000). Here we generalize that approach slightly. Our frame of reference is a deltaic shoreline, but since the model is two-dimensional, it applies to any coastline supplied with sediment via rivers, with no additional transport processes in the offshore region. We propose extensions to the basic model to account for waves and tides, respectively. Figure 2 shows the plan view and cross-section of a two-dimensional delta, at some time t, depositing onto a subsiding basement surface. Under the supply of a unit sediment flux q in (t) at x = 0, this deposit grows with time, its advance measured by the moving position of the shoreline s(t). In developing a model of this planar system, we make only one simplifying assumption: that the foresetthe submerged surface of the delta depositmaintains, through avalanching processes, etc., a constant slope S F .
The starting point in our model development is to consider how the delta shape changes over an incremental time t. In this time increment we see, with reference to Figure 3, that the new area of sediment in the deposit is given by where the prescribed unit volumetric sediment flux at q in (t) includes the contribution of the bed porosity. Figure 3 shows how this new (incremental) area is distributed in the delta. It helps our cause to split the new area into five sub-areas A 1 to A 5 as marked in Figure 3. Area A 1 is the incremental accommodation area created by the relative subsidence under the fluvial topset (above sea level) fraction of the delta at time t: where we denote by σ(x,t) the relative accommodation rate (subsidence plus sea-level change). Area A 2 is the incremental accommodation area created under the submerged foreset fraction of the delta at time t: where we denote by s toe the x position where the constantslope forest intersects the bedrock basement, referred to as the delta toe. Area A 3 is the incremental accommodation area  created under the advanced shorefront. Here it is more difficult to write down the exact area, but we approximate it as where σ rep is a representative value of the rate of subsidence in the vicinity of the delta toe. Note that as our incremental time step becomes small, we expect the contribution from this area to vanish. Area A 4 is the incremental increase in area due to the advance of the shoreline. With our assumption of a constant foreset slope, this can readily be written as where we use D(x = s, t) to represent the depth of the water at the delta toe at time t. With the final incremental A 5 we account for the area of deposition on the fluvial surface: where η(x,t) is the height of the sediment deposit relative to the sea level. Noting that as Δt → 0, ∫ tþΔt t f x; t ð Þdt→Δtf x; t ð Þ, Δs → 0 and Δs Δt → ds dt ,we can substitute our expressions for the incremental areas into Eq. (3) to arrive at the following incremental sediment balance for the delta evolution: To move forward, we recognize that the first integral term on the right represents the partitioning of the incoming sediment to compensate for the creation of accommodation and increase the elevation in the fluvial section of the delta. In this way, we can use Eq. (9) to determine the unit sediment flux arriving at the shoreline x = s: representing the sediment input available to supply the submarine portions of delta and thus advance the shoreline. This definition allows us to manipulate the balance in Eq. (9) in two ways. In the first instance, we can write this equation as which, by the definition of the integral, can also be written as We can always satisfy Eq. (12) by forcing the argument of the integral to be zero at every point in the domain [0, s(t)], that is which we immediately recognize as the Exner equation for the fluvial component of our delta; the appropriate boundary conditions on this equation are We also have the additional condition from the definition of the shoreline flux at x = s: A simple approximation for our problem is to assume that the foreset is a 'cliff face' (i.e. s = s toe ). While obviously physically unreasonable, this greatly simplifies the algebra without changing any of the essential dynamics. In this case our shoreline flux condition becomes which simply tells us that the excess sediment arriving at the shoreline q s advances the shoreline position by filling in the ocean depth. An even simpler additional assumption is to neglect subsidence and assume a horizontal basement. Now our delta is advancing into a water body with a constant depth D and our model for delta growth can be written as subject to the boundary conditions in Eq. (14) and the shoreline conditionobtained from Eq. (15): By making the simple direct substitutions, elevation η for temperature T, sediment flux q for Fourier heat flux q H = À K∂T/∂x and constant depth D for latent heat L, we arrive at the remarkable observation that Eqs (17) and (18) form the classic Stefan problem for tracking the moving phase-change front during the melting of a solid block, a problem that has been investigated intensively in the field of heat transfer for over 100 years. This connection between the flow of heat and of sediment, first made by Swenson et al. (2000), allows us to adapt the extensive techniques and concepts related to melting problems in particular and moving-boundary problems in general to understanding coastal dynamics. It is, in our view, an excellent example of the power of metaphors in science; specifically, we think this metaphor is especially instructive in two ways. First, the metaphor works because of the common mathematics underlying the diffusion of heat to a melt front that advances the front, and the flow of sediment to the shoreline that advances the shoreline. This is despite the fact that the processes are not the samein fact, there is no meaningful sense in which sediment 'diffuses' towards the shore. Rather, the linearized version of gravity-driven fluvial sediment transport leads to an equation that has the same form as a diffusion equation. The second instructive aspect of the melt-shoreline metaphor arises from the shoreline flux condition in Eq. (15), which is a generalization of the so-called Stefan melting condition in Eq. (18). In the Stefan condition, the flux (heat) arriving at the moving front (the water/ice interface) is directly proportional to the speed of advance of the interface; the proportionality constant is given by the latent heat value. Now in real physical systems, the latent heat is effectively a material constanthowever tempted the mathematically inclined might be to explore what could happen if it were not! Adapting the mathematics to the shoreline case via the melt metaphor turns this apparently fanciful idea into reality. In the shoreline condition in Eq. (15), we see that the latent heat termthe depth D(x, t) at the foreset toeis indeed not constant but a function of space and time. In addition we note that, when the foreset has a finite slope, not all of the sediment flux arriving at the shoreline is used towards its advance; a fraction of this flux needs to be used to account for the rate of creation of the relative accommodation space under the foreset. To highlight the fact that the condition in Eq. (15) is a generalized Stefan condition, we refer to it as the shoreline-Stefan condition.

Shoreline Dynamics: Insight from Field and Laboratory
Laboratory studies Our group at the University of Minnesota carried out two explicit experimental tests of the Pitman-Angevine phase-lag prediction for the relation between shoreline and relative sea level (RSL) Paola et al., 2001;Kim et al., 2006aKim et al., , 2006b. In both cases, we compared shoreline response to slow versus rapid change in RSL, and in both cases we were able to constrain all the components of the reference time scale used to nondimensionalize the time scale of RSL change. Thus 'slow' means RSL cycle times that are long compared to the reference time, and 'rapid' cycle times are less than the reference time. It is clear that, at least for the Gilbert-type deltas that can readily be created in the laboratory, the shoreline and RSL change remain nearly in phase for both slow and rapid cycles, and the shoreline response is not strongly attenuated at longer time scales. This is, of course, inconsistent with the prediction of Pitman (1978), though the observations are well described by the simple moving-boundary model summarized earlier. We believe this discrepancy is because the experiments do not satisfy two of the key assumptions in the Pitman formation: (1) that the transport system maintains a surface with constant topographic slope, which is violated in systems with a pronounced foreset; and (2) more importantly, that the sedimentation rate at the shoreline is neither negligible nor constant. The ability of the more general and flexible moving-boundary approach to accommodate such conditions is one of its chief advantages. Below we propose two extensions of the moving-boundary formulation in the hope that these may inspire similar types of experimental testing.

Field studies
To the best of our knowledge there has never been a definitive field test of Pitman's (1978) idea of out-of-phase shoreline motion with respect to changing relative sea level. This is probably due at least in part to the difficulty of knowing the timing of shoreline movement relative to RSL with sufficient precision to measure the phase lag, and of estimating the reference time scale T R needed to calculate the dimensionless cycle time. Given these obstacles to field testing, and that where the model has been tested experimentally its most basic predictions have not been supported, it would be easy to conclude that the Pitman model and its prediction of out-of-phase shoreline behaviour have been relegated to the Museum of Interesting Failed Ideas. In our view, nothing could be further from the truth. The fundamental idea that RSL and shoreline can be out of phase at least implicitly underlies much recent thinking about stratigraphic sequences forced by RSL change. For example, in their comprehensive attempt to synthesize and rationalize terminology across sequence stratigraphy, Catuneanu et al. (2009) indicate (their Figure 15) that the time of maximum flooding (i.e. maximum transgression) coincides with the time of maximum rate of RSL rise, rather than maximum RSL itselfjust as Pitman (1978) originally predicted. They see this as the result of the interplay of rate of change of RSL with sedimentation rate, a different mechanism than that envisioned by Pitman, but still a conceptual cousin to Pitman's proposal. Even more telling is that, upon surveying many years of work relating sequences and relative sea level, Miller et al. (2018) end up with 'we strongly advocate against linking sequence boundaries, MFSs, TSs, and systems tracts to any specific point of relative sea level'. The point is that, while Pitman's model is almost certainly incorrect in some of its basic assumptions and specific predictions, its fundamental contribution of decoupling shoreline motion from RSL change has received the ultimate accolade: to become so embedded in our thinking that we no longer even recognize it.

Extending the Melt Metaphor
Extending the metaphor: waves The above model formulation is restrictive in its geometric treatment of offshore processes and does not allow for the effects of waves on coastal morphology. Can we extend the moving-boundary approach offshore to explicitly incorporate wave-current-driven sediment transport and shallow-marine morphodynamics and, in so doing, expand the metaphor to an offshore system with active transport?
We focus on the problem of how the interplay of river and wave transport controls coastal morphology in the vicinity of river mouths, following the analyses of Wolinsky and Murray (2009) and Nienhuis et al. (2016), but with the aim of linking the shoreline model to the moving-boundary approach to topographic long profiles developed above. To frame the analysis, consider Figure 4, which illustrates in map view coastal-plain, shoreface, and shallow-marine depositional environments on a micro-tidal, constructional coastline under the influence of surface gravity waves. Multiple coastal-plain rivers deliver sediment to a wave-influenced coast, where it is redistributed via a combination of longshore (q ls ) and cross-shelf (q sm ) transport.
A first extension of the basic model might incorporate an explicit treatment of wave current-driven cross-shelf transport but retain the cross-sectional (strike-averaged) formulation. We (Swenson et al., 2005) adopted such an approach in our efforts to model the development of compound clinoforms. We modelled shallow-marine sediment transport as the sum of current-and slope-driven terms [i.e. q sm ¼ f h; ∂h where h is water depth; this formulation generated a nonlinear advectiondiffusion equation governing offshore morphodynamics. The main point relevant to the current discussion is our treatment of the shoreline (shoreface) as a shock condition (i.e. a specified discontinuity in the elevation of the sediment surface): The shock (D) represents physically the depth of the shoreface, which is set by the closure depth for longshore transport. The left-hand side of Eq. (19) is the volumetric rate of sediment stored in the advancing shoreface, whereas the first and second terms on the right-hand side represent, respectively, the fluvial sediment supply to the shoreface and the cross-shelf rate of sediment removal from the shoreface. Justification for treating the shoreface as a shock centres on (1) the short time scale of morphologic adjustment that characterizes this environment and maintains its uniform geometry (e.g. the oft-employed 'Bruun rule'; Lorenzo-Trueba and Ashton, 2014) and (2) the characteristic shoreface slope, which typically is at least an order of magnitude larger than that of the surrounding environments. In effect, we retain a portion of the geometric approach taken earlier, to which we concatenate a second governing equation at larger depths (h > D).
In the context of the melt metaphor, the model developed earlier represents a 'single-phase' Stefan problem, in which one of the two phases (typically the solid) is isothermal and, by extension, conducts no heat. Carrying the metaphor forward, one could interpret the Swenson et al. (2005) model, and the extension we propose here, as 'two-phase' Stefan problems, where cross-shelf sediment transport (q sm ) plays the role of conduction in the second (solid) phase and the shoreface shock (D) represents a constant latent heat.
On many wave-influenced marginsparticularly during times of relative sea-level highstandthe cross-shelf sediment (sand) flux may be relatively small and the marine sediment (sand) budget dominated by longshore transport and construction of extensive strandplain complexes. Such a scenario describes many modern (highstand) coastlines, which are fronted by wide continental shelves largely starved of sediment input. In such systems, especially on sub-Milankovitch time scales, it is difficult to ignore the three-dimensionality of sedimentation. We thus propose a simple extension of the Swenson et al. (2005) model in which strandplain development is coupled explicitly to the feeder coastal-plain rivers on time scales where subsidence is unimportant. We borrow from well-established, 'single-equation' models of both fluvial and shoreface morphodynamics of comparable sophistication. The emphasis here is less on detailed treatments of fluvial and longshore sediment transport and more on the careful coupling of the depositional environments at the channel termini, which we treat as moving boundaries.
In map view, we denote by s(y, t) the lateral position of the shoreline (shoreface) with respect to the y axis ( Figure 4). Assuming the representative morphology-shaping wave field is described by low-angle waves, then, to a good approximation, the longshore sand flux (q ls ) varies linearly with the plan-view gradient of the shoreline: where κ is a diffusivity that varies nonlinearly with the representative wave height (Komar, 1983;Ashton et al., 2001). The shoreline evolves in response to divergence in longshore sand flux: Equation (21) holds everywhere but the distributary channel termini, which are located at positions y i . The shoreline is thus piecewise continuous, with jump discontinuities in shoreline gradient (our proxy for longshore sand flux) at the termini (y i ) of the feeder channels.
We treat fluvial morphodynamics as before: the long-profile evolution of the i th coastal-plain river obeys a transient linear diffusion equation: where η i is the long-profile elevation, x i is a streamwise coordinate and υ i a diffusivity that varies linearly with water discharge (Paola, 2000). Boundary conditions on Eq. (22) are a specified sediment flux (q soi ) at x i = 0 and a specified elevation (sea level; η i = 0) at the channel terminus (x = s ci ). Longshore transport couples to fluvial input at the channel termini. For simplicity, we assume the channels are situated in belts of representative width W i , where the W i are small relative to the spacing of channels. The coupling condition equates the rate of channel progradation to the difference between fluvial sediment supply reaching the channel terminus and extraction by longshore transport: where D again is the closure depth for longshore transport. The left-hand side of Eq. (23) represents sediment sequestered in the prograding channel belt; the first term on the right-hand side is the fluvial (volumetric) sediment discharge at the terminus; and the third term captures wave-driven extraction of sediment on both sides (À, +) of the channel. Equation (23) serves as an additional Neumann boundary condition on both Eqs (21) and (22), thereby leaving both governing equations overspecified and, in turn, allowing for determination of the seaward position of the channel terminus.
As an application of the above theory, consider a hypothetical highstand coastline receiving input from a pair of coastal- plain rivers ( Figure 5). The sediment supplies (Q s ) to the left and right rivers are 2 × 10 6 and 10 6 m 3 a À1 , respectively; the fluvial (υ) and coastal (κ) diffusivities are 8 × 10 6 and 2 × 10 5 m 3 a À1 ; both rivers possess a channel-belt width (W) of 500 m; and the closure depth for alongshore transport is 5 m. Figure 5 shows model predictions for one millennium of evolution. The right river remains active throughout the entire model period, whereas the left riverwith twice the sediment dischargeis abandoned at t = 500 a. In map view, alongshore transport of sand from both rivers constructs characteristic cuspate strandplains (Figure 5a), with the left river-associated cusp being reworked extensively following abandonment.
A key aspect of this model is its rigorous coupling of channel termini to adjacent shorefaces via Eq. (23), which is applied at each channel terminus. Figure 5b shows the trajectories of the left and right-river termini throughout the millennium of coastal evolution; Figures 5c and d show the corresponding fluvial long profiles. In this strongly wave-influenced system, strandplain construction is the dominant sink for sand and, consequently, progradation of the channel belts is sufficiently slow to allow the fluvial long profiles to maintain a quasi-equilibrium ('graded') state, which, by Eq. (23), is linear. The trajectories in Figure 5b reflect the (1) interplay of sediment deposition in the prograding shoreface and aggrading and prograding channel belts and (2) interferencevia longshore transportof fluvial systems.
Scaling the above governing equations and closure condition [Eqs (21)-(23)] yields a single, dimensionless number ( D ffiffiffiffiffi κυ p =Q so ) that controls the extent of coupling between channel progradation and longshore transport. When D ffiffiffiffiffi κυ p =Q so < 1, sediment partitioning between channels and shoreface is comparable in magnitude. In contrast, if D ffiffiffiffiffi κυ p =Q so >> 1, then the bulk of sand supply is sequestered in beach ridges, and the feeder distributary channels are slaved to longshore transport. The system in Figure 5 is characterized by D ffiffiffiffiffi κυ p =Q so~6 .
While strandplain deposition dominates the budget here, the importance of careful coupling to the feeder fluvial systems should not be downplayed, as it is the contemporaneous responses of these environments to sea-level change that generate key stratigraphic surfaces (e.g. wave ravinement surfaces and incised valleys). In that regard it would be useful to explore ways to incorporate into this model insight from more detailed models of river mouth-wave interaction, such as the recent work of Nienhuis et al. (2016).
Stronger coupling between channel and shoreface dynamics is more likely at smaller scales (e.g. where distributary channel density is greater), such as a young (sensu Jerolmack and Swenson, 2007) delta developing in an environment with moderate wave energy. The Tinajones delta of northwestern Colombia provides one such example. Here, in approximately 1943, the Sinu River avulsed, spanned a narrow beach to the Caribbean Sea and nucleated a new delta (Troll and Schmidt-Kraepelin, 1965), which then underwent a systematic evolution towards a more wave-influenced state (Figure 6 and imagery in Suarez, 2004). In the first few decades of development (Figure 6a), the foreset area of the Tinajones delta was relatively small, ensuring that fluvial delivery was sufficiently large to maintain strong coupling between distributary channels and wave-driven sediment transport (D ffiffiffiffiffi κυ p =Q so < 1 in the model above). However, with increasing topset/foreset area, and an additional distributary channel, wave-driven transport began to dominate, and the Tinajones delta was driven to a state of stronger wave influence (D ffiffiffiffiffi κυ p =Q so >> 1), as reflected in its overall appearance (Figure 6d). We applied the above theory to a simplified version of this scenario, in which a pair of closely spaced distributary channels receive equal sand discharges. We modelled 20 years of delta growth for two scenarios. In the first, meant to represent an 'early' phase of delta growth, we employed a large sediment discharge (Q so = 10 6 m 3 year À1 ) and weak wave climate (κ = 10 4 m 2 year À1 ) such that D ffiffiffiffiffi κυ p =Q so = 0.5; in the second, which reflects a more mature delta, we halved the sediment supply (Q so = 5 × 10 5 m 3 year À1 ) and increased the wave strength to generate an order-of-magnitude greater wave influence (i.e. D ffiffiffiffiffi κυ p =Q so = 5). We held other variables constant: υ = 10 4 m 2 a À1 , W = 200 m and D = 5 m. Figure 6 summarizes the results of this modelling experiment. The 'early-phase' delta with minimal wave influence partitions significantly more of its sediment supply to channel-belt progradation/aggradation (Figure 6e), resulting in significantly greater plan-view coastline rugosity associated with the seaward extension of distributary channels (Figure 6b). Continued growth leads to relative strengthening of the wave influence and hence a smoother shoreline (Figure 6f). This is also visible in the Tinajones field example (Figure 6d), where the individual distributary tips display less pronounced protrusion from the surrounding coastline than did the newly avulsed channel (Figure 6a). Also, in contrast to the fluvial long profiles in Figure 5, which reflect a state of quasi-graded behaviour, the long profiles of the 'early' system here (Figure 6c; D ffiffiffiffiffi κυ p =Q so = 0.5) show considerable curvature that is indicative of nongraded (disequilibrium) sedimentation.
Extending the metaphor: tides and mushy boundaries One of our main themes in this review has been using the analogy of a Stefan melting problem to model and understand shoreline dynamics. How far can we push this analogy? The Stefan melting problem is a particular example of a phase-change problem. The key feature in the problem is that melting (the phase change from the solid to the liquid Figure 6. (a) Tinajones delta, Colombia, in 1957, approximately 15 years after its nucleation; line drawing modified from Troll and Schmidt-Kraepelin (1965). (b, c) Modelled evolution of shoreface and channel mouths (map view) and fluvial profiles (cross-section) over a 20-year interval in a deltaic system with strongly coupled (D ffiffiffiffiffi κυ p =Q so = 0.5) shoreface and channel dynamics; note curvature of fluvial profiles. (d) Tinajones delta in 1989; note morphological features indicative of increasing wave influence, relative to 1957; air photo modified from Suarez (2004). (e) Comparison of channel-terminus progradation in strongly and weakly (dashed) coupled systems. (f) Evolution of shoreface and channel mouths in a less strongly coupled deltaic system (D ffiffiffiffiffi κυ p =Q so = 5); compare to (b). state) occurs at a unique temperature. However, in more general phase-change systems (e.g. with impurities), the phase change can occur across a range of temperaturesa socalled 'mushy phase change'. Can we find an analogy for the mushy phase change in delta shoreline buildinga mushy shoreline? Remembering that in the melt metaphor the temperature is represented by surface elevation η, we see that a mushy shoreline is, in effect, one in which the transition from land to water occurs over a range of elevations, rather than a single one, and hence over a range of x values: a 'shore zone' rather than a shoreline. This rather nicely describes a shoreline influenced by tides, for which sea level (i.e. the melting temperature) varies over a delimited, small range, on a time scale (about half a day) that is much shorter than that of the morphodynamic evolution we are interested in tracking.
To apply the mushy-boundary idea to a tidal coast, let us consider a long-profile cross-section building into an ocean with a fixed depth L (Figure 7). Let us assume that above a particular subaerial elevation η * > 0 on this profile the supply of sediment is used exclusively to build the land deposit. This elevation thus marks the high-tide line. Below this elevation, however, we will assume that, in plan view, the surface is made up of both land and water areas ( Figure 7). As we move seaward, the water area (averaged in the lateral direction) increases from a value of zero to the point, at and beyond η = 0 (the low-tide level), where there is only water in plan view. The elevation range η * ≥ η ≥ 0 defines the intertidal zone, and hence our mushy shoreline, within which the varying fraction of land vs. water can be thought of as a bulk representation of the network of tidal channels. We assume that any sediment entering this region is partitioned into two fractions. One fraction builds the subaerial elevation; the other fills in the water depth. We can track this sediment by defining a total deposited sediment thickness as where L À L 0 > L m > 0 is the average depth of the submarine deposit at a point within the intertidal zone, L(x,t) is the depth from the sea level to the basement and L o (x,t) is the difference (drop-off) between the tidal-flat deposit and the basement at the extremity of the intertidal region. For a simple illustrative analysis, we assume that the basement and drop-off depths are constant and that the deposit depth in the intertidal zone can be modelled as Taking the sediment flux as proportional to the fluvial slope allows us to write down the following model for tracking the growth of the delta: where q in is the unit discharge of sediment introduced at the origin and the relationship between H and η is given by otherwise 2 6 6 6 6 6 6 6 6 6 6 6 4 The value ν is a diffusion coefficient: at elevations above η * , we assume that this value is a constant ν up while, to model the effect of tides, in η * > η > 0 we will set ν ¼ increase in diffusivity embodied in Eq. (27) represents the effect of adding ebb-tidal discharge to the input discharge in moving sediment; it has the effect of decreasing the slope in the tidally influenced region as observed in, for example, the tidally influenced region of the Ganges-Brahmaputra-Meghna Delta of Bangladesh (Wilson and Goodbred, 2015).
Returning to the melt metaphor, we note that this model is essentially the well-known enthalpy method (Voller and Cross, 1981;Crank, 1984) which has previously been applied to the case of sharp shorelines (Voller et al., 2006;Lorenzo-Trueba and Voller, 2010). The real convenience of the model is that it is applicable throughout the entire domain x ≥ 0 and can be solved using a simple explicit-time-intergration finite-difference scheme.
An example schematic calculation of a profile (q in = 1,L = 1, L m = 0.9,L o = 0.1,η * = 0.2,ν up = 50,ν down = 200) is shown in Figure 8. This profile represents a lateral average over topography. Figure 8 shows the predicted subaerial surface η in the tidal region (e.g. the tops of tidal marshes) and also an average surface Figure 7. Definition sketch for a simple mushy-boundary model for a tidal coastline. elevation for the sediment deposit H À L that accounts for both marsh tops and tidal-channel bottoms.
At this point we stress that this model is a first step that illustrates how the analogy of the Stefan problem of melting can be extended to model tide-influenced shorelines. We believe this approach can eventually be merged into a general, integrated 1-2D treatment of river, tidal and wave controls of shorelines.

Metaphorical Shorelines
As important as the shoreline is, it is hardly the only major facies boundary one encounters in the source-sink system. The moving-boundary approach we advocate to model it can be applied to any dynamic facies or environmental boundary separating different transport regimes. For example, Figure 9 shows two examples of multiple interacting alluvial systems, one experimental and one from the field. Both images show river systems with multiple sources of water and sediment: one throughgoing axial channel and a series of subordinate transverse systems that feed sediment into the axial system. Such systems are well known from rift fills such as the Rio Grande Rift of the southeastern USA (Mack and James, 1993;Leeder et al., 1996). The first key point for our purposes is that despite their familiarity across the Earth's surface, there is nothing 'natural' or obvious about this arrangement of the fluvial system; rather, the geometry results from self-organization of the variously sourced systems interacting with each other and the subsidence field, mediated by a set of dynamic moving boundaries (Grimaud et al., 2017). A second key point is that from the standpoint of the transverse systems, the axial system represents a substantially different transport environment, typically with larger transport capacity associated with higher discharge. The analogy between the axial-transverse boundary and shorelines is particularly clear, a point first recognized by Kim et al. (2011). The entrainment of transversely sourced sediment into the axial system ('toe cutting' of Leeder and Mack, 2001), which exerts an important control on the boundary, is formally analogous to the entrainment of river sediment into the longshore transport system at the coast, as discussed above. For example, as detailed in Kim et al. (2011), when the transport capacity of the longitudinal system increases, it moves the transverse-longitudinal boundary up into the transverse fan. This is in effect a transgression, analogous to shoreline transgression caused by wave erosion as, for example, that occurring now on the Canterbury Plains, South Island, New Zealand (Leckie, 1994). They may be far from any lake or ocean, but such 'shorelines' are dynamic boundaries that migrate in time and control the partitioning between axial and transverse facies in the deposit. Though it has even less resemblance to an actual shoreline, the start of deposition in a river systemthe basement-alluvial transitionis likewise an abrupt change in transport regime (Kim and Muto, 2007), and thus amenable to a shoreline-like analytical treatment (Lorenzo-Trueba et al., 2009). The most striking finding to emerge from the work of Kim and Muto (2007) is that in some ways this 'shoreline' behaves inversely to the actual shoreline: where the latter shows stronger autogenic fluctuations when being driven against the sediment transport direction (i.e. transgression and rising RSL), the basement-alluvial boundary fluctuates most during regression of the overall system (Kim et al., 2006a, Kim et al., 2014. Another abrupt change in transport regime in river systems occurs at the gravel-sand transition in alluvial rivers (Ferguson, 2003;Venditti and Church, 2014;Lamb and Venditti, 2016;Delorme et al., 2017). In this case the regime change is from nearthreshold conditions in the gravel case to transport well above threshold with a significant fraction of the bed-material load in suspension in the sand-bed case. This was formulated in moving-boundary form by Marr et al. (2000). The behaviour of this boundary is more straightforward than that of the bedrock-alluvial case, in large part because it is controlled directly by mass balance: as long as abrasion does not play a major role in limiting gravel transport (Paola, 1988;Parker, 1991aParker, , 1991bDingle et al., 2017), it occurs when the supply of gravel has been exhausted by deposition. In many rivers this occurs over a relatively small fraction of the overall system length (Ferguson, 2003), so to a first approximation the moving boundary is indeed sharp. However, the 'mushy-boundary' method described above can readily be adapted to account for the transition region. Using the modelling techniques in Marr et al. (2000), let us consider, for illustrative purposes, building a fluvial fan with a planar cross-section onto a nonsubsiding horizontal basin aligned in the x direction. The sediment flux into the system, q in (sediment volume + porosity) is a mixture of gravel and sand, with a gravel volume fraction of 0 < α g < 1. For this case, a simple but reasonable model for the sediment flux of the gravel or sand fractions at a point on the fan surface is where the subscript f = g or s indicates the sediment component (gravel or sand), ν f is a diffusivity, S f is the threshold slope for transport for a given component, and η(x,t) is the current sediment height (gravel + sand) of the fan. The flux model in Eq.
(28) essentially states that a sediment component will be transported at a given point only when the slope of the fluvial surface is greater than a threshold value and that this flux cannot be larger than the fraction of this component in the input flux. With this constitutive relationship, the building of the fan can be modelled through a numerical solution of the Exner balance, in this case equating, at each point in the domain, the transient increase in the sediment height of a component, η f , to the divergence of its flux. Note, for convenience of computation, that we set equal diffusivities and use an arbitrary choice of space and time dimensions; the phenomenological drivers in this problem are the threshold slopes and gravel fraction.
As an example of what output might be expected from such a model, we consider a situation where the threshold slope for gravel is larger than the threshold slope for sand, S g = 0.05 > S s = 0.025. In this case we would expect the gravel to be deposited in the upstream portion of the fan and the sand in the downstream. At the input (x = 0), where the fan slope is steepest, the sand will bypass without depositing. Thus, the gravel flux component sets the initial slope on the fan. Here we set diffusivity and input flux so that the initial slope of the fan is~1.7 times the gravel threshold.
At any point in time there are two points of interest on the fan surface, the point where the gravel in the feed runs out (i.e. the point where all the gravel has been deposited) and the point where the sand begins to deposit (i.e. the point where the sand flux calculated for the current slope is less than the sand flux in the input). The relative position of these two points along the fan surface determines if the gravel-sand transition is sharp or 'mushy'.
The left panel in Figure 10 shows the stratigraphy of the fan, for a given point in time, when the fractions of the input flux are α g = 0.75,α s = 1 À α g = 0.25, respectively. In this case, the gravel runs out before any sand deposits. This results in a sharp transition between the gravel and sand domains in the fan stratigraphy and a slope break in the surface profile at the transition; in order to deposit sand after the gravel runs out, a shallower slope is required. As the sand fraction in the feed is decreased further, the transition remains sharp and the slope break becomes more pronounced. By contrast, when the feed sand fraction is increased, the slope break at the sharp transition diminishes. In fact, there will be a point where the slope is continuous across the transition. We can see, by matching the slope at the gravel run-out (S g ) with the slope at the first sand deposition ( q s νs þ S s ), that this slope continuity occurs when the sand fraction reaches the value which for our current example is α s = 0.375. When the sand fraction exceeds this value, we would expect the sand deposition point to be above the gravel run-out point. Between these points we establish a 'mushy region' in the stratigraphy, a region where both gravel and sand are simultaneously deposited. The right panel in Figure 10, using the same input flux but gravel and sand fractions α g = 0.5,α s = 0.5, shows the appearance of this mushy region. In this case the fan surface is continuous through the transition, with no slope break. While these are early results, they demonstrate the possibility of evolving mushy-region transitions in landscape systems. In particular, the criterion in Eq. ((29)) provides a testable hypothesis for the appearance of mushy gravel sand transitions.
Stepping back, it is not hard to see the sequence of major environmental transitions going from source to sink down the sediment-routing system (Allen, 2017) as a sequence of moving boundaries, amenable to treatment in a common framework, with differing linking conditions at each boundary reflecting the handoff dynamics at the boundary between transport laws within the two communicating domains. The essence of the common framework is the structure outlined earlier, with examples of how the framework can be adapted to specific transport conditions and boundary characteristics (wave and tidal coastlines, and the fluvial gravel-sand transition). But the larger message is that in a sense, all major environmental facies boundaries can be seen as, in essence, shorelines.

Conclusions
1. Pitman's (1978) demonstration that the apparently obvious link between shoreline and sea level, based on common experience, can break down on planetary time scales was a conceptual breakthrough that underlies shoreline modelling and sequence stratigraphy to this day. 2. The metaphor between the melt front between liquid and solid phases in a heated solid and the shoreline of a fluvial coastal plain provides the basis for a rich theory of shoreline dynamics. 3. Extension of the melt metaphor to waves leads to a 'twophase' shoreline model in which wave action results in local transgression once channels are abandoned and smoothing as a new delta lobe grows into a basin with wave transport. 4. Melting of alloy materials in which the phase change occurs over a range of temperatures yields a so-called 'mushy' melt front that provides the framework for a 1D model for a shoreline with an intertidal zone, and a model of a gradual gravel-sand transition in rivers. 5. The shoreline is the type example of a geomorphic interface that can be successfully treated using moving-boundary methods, but the technique can be applied to environmental interfaces across the entire source-sink system: a set of environments linked by metaphorical shorelines. in this case the gravel (red) runs out before the sand (yellow) deposits, giving a sharp transition and a slope break in the fan surface. Right: results for sand fraction = 0.5; now there is a 'mushy' region in the cross-section (white space), where both gravel and sand are deposited together.