A Time-Continuous Embedding Method for Scalar Hyperbolic Conservation Laws on Manifolds

A time-continuous (tc-)embedding method is first proposed for solving nonlinear scalar hyperbolic conservation laws with discontinuous solutions (shocks and rarefaction waves) on codimension 1, connected, smooth, and closed manifolds (surface PDEs or SPDEs in R2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbb {R}}^2$$\end{document} and R3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbb {R}}^3$$\end{document}). The new embedding method improves upon the classical closest point (cp-)embedding method, which requires re-establishments of the constant-along-normal (CAN-)property of the extension function at every time step, in terms of accuracy and efficiency, by incorporating the CAN-property analytically and explicitly in the embedding equation. The tc-embedding SPDEs are solved by the second-order nonlinear central finite volume scheme with a nonlinear minmod slope limiter in space, and the third-order total variation diminished Runge-Kutta scheme in time. An adaptive nonlinear essentially non-oscillatory polynomial interpolation is used to obtain the solution values at the ghost cells. Numerical results in solving the linear wave equation and the Burgers’ equation show that the proposed tc-embedding method has better accuracy, improved resolution, and reduced CPU times than the classical cp-embedding method. The Burgers’ equation, the traffic flow problem, and the Buckley-Leverett equation are solved to demonstrate the robust performance of the tc-embedding method in resolving fine-scale structures efficiently even in the presence of a shock and the essentially non-oscillatory capturing of shocks and rarefaction waves on simple and complex shaped one-dimensional manifolds. Burgers’ equation is also solved on the two-dimensional torus-shaped and spherical-shaped manifolds.


Introduction
Partial differential equations on a surface (surface PDEs or SPDEs) have received growing interest over the last years. Problems of interest areas arise in a wide range of fields such as image processing (e.g., the restoration of a damaged pattern on a surface [4] and image of the human brain [24]), biological systems (e.g., study domain formation in vesicles [2] and wound healing process [27]), mathematical physics (e.g., flows and solidification on surfaces [25,26]) and computer graphics (e.g., real-time fluid visualization on surfaces [1]).
Instead of studying linear and quasi-linear SPDEs with smooth solutions ( [12,14] for level set approach, [31] for a semi-Lagrangian method, [37] for meshfree finite difference method, [3,7,10,22,[28][29][30] for closest point (cp-)embedding methods) on Cartesian grids and orthogonal gradient grids, as a proof of concept, this study explores the application of the embedding methods and develops a numerical scheme, based on our extensive knowledge of essentially non-oscillatory shock-capturing schemes, to solve the (nonlinear scalar) hyperbolic conservation laws (Burgers' equation) with a discontinuous solution on one-dimensional manifolds. For the rest of the paper, the manifold is implicitly referred to as the one-dimensional, connected, smooth, closed, and complex (or simple) manifold unless specified otherwise.
After reviewing some necessary definitions and assumptions in Sec. 2, we propose a new time-continuous (tc-)embedding method based on the classical closest point (cp-)embedding in [30] in Sect. 3. The new tc-embedding method improves upon the classical cp-embedding method, which requires re-establishments of the constant-along-normal (CAN-)property of the extension function at every time step, in terms of accuracy and efficiency by incorporating the CAN-property analytically and explicitly in the embedding equation. The tc-embedding SPDEs are solved by the second-order nonlinear central finite volume scheme (CFV) with a nonlinear minmod slope limiter in space, and the third-order total variation diminished Runge-Kutta scheme (RK-TVD) in time. An adaptive nonlinear essentially non-oscillatory polynomial (ENO) interpolation is used to obtain the solution values at the ghost cells required by the CFV scheme. The main contributions of this study are in the following two areas. They are i) a novel time-continuous (tc-)embedding that holds for all time continuously in contrast to the classical cp-embedding that is only true at a fixed time (see Theorem 1 for the formal statement), and ii) the first attempt in solving hyperbolic SPDEs with discontinuous solutions on a manifold. In Sect. 4, we detail a numerical recipe (CFV, RK-TVD, and ENO interpolation) for solving the resulting embedded SPDEs. In Sect. 5.1, accuracy and efficiency tests of the tc-embedding method versus the cp-embedding method are given. The results indicate that the tc-embedding method has a lesser numerical dissipation, resolves the fine-scale structures more efficiently, and is roughly two and half times faster CPU times than the cp-embedding method. In Sect. 5.2, the Burgers' equation with convex flux, the traffic flow problem with concave flux, and the Buckley-Leverett equation with non-convex flux demonstrate the essentially non-oscillatory shock-capturing property and the efficient resolution of fine-scale structures of the tc-embedding method on a complex manifold. In Sect. 5.3, numerical results are also given for extending the tc-embedding method for solving the Burgers' equation on two-dimensional torus-shaped and spherical-shaped manifolds embedded in R 3 . Conclusion and future works are given in Sect. 6.

PDEs on Manifolds
Let M ⊆ R 2 be a (one-dimensional, connected, smooth, and closed) manifold with welldefined unit outward normal vector and tangent space, denoted byn ∈ R 2 and T M , respectively, for every point p ∈ M.
Also, let Q M : M × R + → R be some conserved variables associated with a smooth flux function F : R → R. In this paper, we consider SPDEs of scalar hyperbolic type for the unknown conserved function Q M ∈ H 1 (H 2 (M); [0, T ]) in the form of for some known vector field v M : M → T M . The surface gradient operator ∇ M in (1) can be defined as in terms of the standard Euclidean gradient ∇ in R 2 via projections to T M .

Embedded PDEs
For some δ > 0, we begin by defining the narrow band domain around the manifold by The goal of this section is to recast the well-posed PDEs on M in (1) in . Firstly, we need the relationship between surface and Euclidean gradient operators on time-independent functions on the manifold and their extensions.
Proposition 1 gives a strong motivation to extend manifold functions Q M to Q E in such a way that the normal-derivative ∂n Q E remains zero everywhere on M. Note that this is not yet a unique extension operator.
Before proceeding, let us define the closest point mapping cp : on in (3) with some sufficiently small δ [30]. The existence of such , in which cp is well-defined, is guaranteed provided that M has enough smoothness; see [8]. For any (time-independent) manifold function Q M ∈ M → R, the constant-along-normal (CAN) extension is given by whose function value is constant along with normal directions of M.
Although the notation difference between (7) and (12) is small, the implication is huge in the computational aspect. Instead of embedding the original problem (1) in the form of (7) at every time step, we only need to tc-embed (1) once and solve (12) until we reach the final time T . It is because, as shown later, the extended velocity field v E depends exclusively on the geometric properties of the manifold, which is time-independent in this study. One can easily see that the requirement (9) implies The realization of the proposed tc-embedding hence lies on finding an appropriate extended velocity field v E : → R 2 to the vector field v M : M → R 2 . To do so, we define the signed-distance function dist : Now, the following theorem provides the mean to adjust the extended velocity field v E so that the CAN-property of Q E is preserved over time. (3)

Theorem 1 Suppose that M is smooth enough and in
where κ(cp(x)) andκ(cp(x)) denoting the absolute and signed curvature of M at cp(x), respectively, finite and exist, then Proof Due to the CAN-property, it suffices to show that for any fixed x ∈ ⊃ M. Note that the right-hand side of (14) is a v M -directional derivative of q at cp(x) ∈ M. Consider a (local) arc-length parameterization γ (s) of the manifold M defined by its arc-length variable s ∈ N (0) in some -neighborhood around the origin with γ (0) = cp(x) and orientation given by the unit tangent vector d ds Obviously, v E must point to theT direction. We only need to study the scaling effect when the point of evaluation in the directional-derivative ∂ ∂T q is shifted from cp(x) to x.
Let D := dist x, cp(x) , that is a fixed real number because x is considered as fixed. Since the outward unit normal vectorn(cp(x)) =n(γ (0)) of M depends on the orientation of the manifold (definition of counterclockwise) and the normal vectorN(0) obtained from γ (s) depends on the orientation of the curve γ (s), we have We begin by considering the directional derivative as a derivative on γ (s), i.e., Gateaux derivative, and shifted the point of evaluation using the CAN-property as follows: The last equality uses another definition for directional derivative based on the rate of change in the directionT ± DN , i.e., Fréchet derivative. Using the Frenet-Serret formula and the fact that M is planar, we find that The scaling factor we desired is given by (16) at s = 0, and hence, the formula for v E in (13) is deduced.
In the proof above, we have used the definition of the absolute and signed curvatures, that is, (13) is time-independent. It can be pre-computed once at each point x in the initialization stage of the algorithm and used subsequently in the time-stepping scheme. Now that the embedding process is completed solving the SPDEs (1) can be transformed into solving the tc-embedded PDEs (12) in a narrow computational domain with any suitable spatial solver in the two-dimensional Euclidean space. In the coming sections, we will describe the tc-embedding method for solving (12) based on the governing equations and the CAN-property in its solution, i.e., instead of the zero Neumann boundary conditions. Using this alternative approach, we propose an easy recipe to propagate non-smooth solutions on one-dimensional manifolds.

Implementation of the Tc-Embedding Method
This section will describe the numerical method in solving tc-embedded PDEs (12) with the spatial operator L (11) in a computational tube of radius δ.
We shall denote the computational domain in the two-dimensional Cartesian domain as (bounded within two pink lines) where the one-dimensional manifold M (black line) is embedded inside as illustrated in Fig. 1 with a circular manifold. With the Cartesian domain , we discretize the domain with uniformly spaced computational cells with cell sizes x and y in the horizontal x-direction and the vertical y-direction, respectively. In Fig. 1, the cell centers x = (x i , y j ) are shown. Following the suggestion in [28], we restrict the computational tube radius δ to be where the critical radius δ * satisfies the condition δ * max{κ(x) : x ∈ M} < 1, where κ is the curvature of M. Since we will be using a second-order finite volume scheme (CFV) and Fig. 1 Simple illustration of the distribution of data points X = X I X G ⊂ in the computational tube with the sets of interior X I (blue) and ghost X G (red) points. The black line represents the manifold M (Color figure online) a second-order ENO interpolation, to solve (12) at a given x = (x i , y j ) (excluding the ghost cells), the computational stencil is p = 2 cells in all directions. This lower bound of δ ensures that there is always a sufficient number of computational cells available in the neighborhood of a given cell center x to perform the necessary spatial operations given the order of the numerical scheme p.
Together, as illustrated in Fig. 1, we define a set of uniform spaced cell center (grid) points (x i , y j ) with cell size x and y in by contained inside the computational tube of radius δ. The set X is further sub-divided into two subsets; • The set of interior cells (blue solid symbol) X I := {x ∈ X : for computing the flux gradient in(11)}, • The set of ghost cells (red hollow symbol) X G := X \ X I .
Ghost cells X G are near-boundary cells that cannot be updated due to missing stencil's information, but are needed to compute the flux gradients in (11) for those cells at the edges of X I . The size of X G generally grows linearly as the spatial order of the scheme. For the second-order scheme used here, one cell is needed inside the domain boundary. At any given time t ≥ 0, we can use the available numerical solution Q E (X, t) at X to compute the numerical flux gradient at X I by which will be discussed in greater detail below. Together with the extended velocity field v E = (u, v) defined in (13), the semi-discrete approximation L N of the spatial operator L defined in (11) becomes Finally, we obtain the semi-discrete system of ODEs where the initial condition is specified by applying the cp-mapping to the initial conditions in (1). We note that while the system of ODEs is evaluated only at the interior points X I , the initial condition must be made available in the whole computational tube X that includes the ghost cells X G . The numerical solution Q E (X I , t) of (20) can be solved by advancing in time via the explicit third-order total variation diminishing (TVD) low-storage Runge-Kutta (RK-TVD) scheme [13] from time t = t n with the discrete solution Q n I = Q n E (X I ) = Q E (X I , t n ), to time t = t n+1 = t n + t with a time step t, Q n+1 with the given initial condition Q 0 are the intermediate solutions and times (t 1 = t n + t, t 2 = t n + t/2) at the kth RK-TVD stage. The time step t at each RK-TVD step is determined by the CFL stability condition |∂ F/∂ Q| is the maximum characteristic speed of the hyperbolic conservation laws on X. The CFL stability condition for an explicit numerical SPDE solver depends not only on the time step t, grid spacing r , and spectral radius of the Jacobian of the flux, but also on the geometric properties, namely, curvature κ expressed via v E (13), of the given manifold. For a manifold with a large curvature κ, the time step t will be reduced accordingly. In the following discussion, for simplicity, we shall take x = y and CFL = 0.45 unless stated otherwise. Readers are reminded that t k , k = {1, 2} are the intermediate times between the times t n (at the time step n) and t n+1 (at the time step n + 1).
Since the tc-embedded PDEs (9) with the extended velocity field v E satisfies the CANproperty intrinsically, the Neumann boundary condition in (12) is discarded in exchange for finding the solution Q E (x, t) at the ghost cells x ∈ X G .
Given the numerical solutions {Q , the two missing components of the numerical scheme in the discussion are: 1. At the beginning of each RK-TVD stage with k = {n, 1, 2} and for all x ∈ X I , computing the flux gradient ∇ F(Q E (X I , t k )) using the shock-capturing central finite volume (CFV) scheme with the nonlinear minmod slope limiter (see Sect. 4.1), and 2. At the end of each RK-TVD stage with k = {1, 2, n + 1} and for all x g ∈ X G , updating the numerical solution for the ghost cells Q E (X G , t k ) with the essentially nonoscillatory (ENO) interpolation at ξ = cp(x g ) ∈ M, which satisfies the CAN-property, based on a few specially selected Q E (X I , t k ) in the closest neighboring cells around ξ (see Sect. 4.2).

Nonlinear Construction of the Flux Gradient ∇F(Q E (X I , t k ))
At any given time t = t k , k = {n, 1, 2}, the numerical flux gradient , at a given cell center x = (x i , y j ) ∈ X I is approximated via a second-order shock-capturing nonlinear central finite volume (CFV) scheme with the nonlinear minmod slope limiter [19] in the two-dimensional Euclidean space, which is a class of nonlinear numerical schemes designed for solving nonlinear hyperbolic partial differential equations. Since the computational stencils X are embedded in the Cartesian rectangular domain with uniform cell-sizes x and y in both x-and y-directions, respectively, the numerical approximation of the flux gradients can be carried out simply in a dimension-by-dimension order. Without loss of generality, we will only briefly describe a generic numerical framework in computing the x-component ∂ F ∂ x via a one-dimensional CFV scheme by assuming the solution Q E and the flux F are functions of x alone and independent of y. The same numerical framework for computing the y-component ∂ F ∂ y can be carried out similarly by assuming the solution Q E and the flux F are functions of y alone and independent of x. Therefore, we will describe the one-dimensional CFV scheme's generic numerical framework without the burden of irrelevant superscript and subscript notations needed in a multi-dimensional case. Furthermore, since we are only interested in finding the spatial gradient of the flux at a given time t, the time-dependency of the solution Q E and the flux F and their corresponding subscript and superscript notations are also ignored in the discussion below.
In a one-dimensional generic numerical framework of the CFV scheme, at a given cell ] with the cell center x i , cell boundaries x i± 1 2 , and cell-size x, the computational stencil {I i−1 , I i , I i+1 } and their corresponding numerical solution {Q i−1 , Q i , Q i+1 } are employed to compute the gradient of the flux written in the conservative form at x i exactly as to guarantee the numerical solution of the CFV scheme is a weak solution of the hyperbolic conservation laws. The numerical flux F i+ 1 2 at the cell boundary x i+ 1 2 can be approximated by the global Lax-Friedrichs (LF) flux as where the slopes are approximated as Due to the non-linearity of the hyperbolic PDEs, a finite-time singularity (discontinuity) might formed. To mitigate the Gibbs oscillations from polluting the solution and destabilizing the numerical scheme, the slope (Q x ) i should be limited. Therefore, we apply the generalized nonlinear minmod slope limiter [19] ( where the general multivariate minmod function is In this study, unless otherwise specified, we take ω = 1.25 for all the numerical examples.

Essentially Non-Oscillatory (ENO) Interpolation for Ghost Cells X G
At any given time t = t k , k = {1, 2, n + 1}, the numerical solution at all ghost cells Q E (X G , t k ) are needed to be specified to allow an update of the numerical flux gradient in the interior cells ∇ F(X I , ·) in the next RK-TVD stage/step. Using the CAN-property of the embedding method and the closest point mapping cp(x), the value at the end of the RK-TVD stage/step. Without loss of generality, we shall ignore the time-dependency of the variables and their corresponding superscript notations in the discussion for simplicity. One can argue that a linear (e.g. central polynomial) interpolation method should work well if the underlining function is sufficiently smooth. However, one cannot assume any smoothness of the nonlinear hyperbolic PDEs' solution in general. The solution can become discontinuous at any spatial location and time. Hence, one should not use the linear high order interpolation method near a discontinuity lest the interpolation across the discontinuity will result in numerical oscillations. To get around this predicament, a nonlinear interpolation method that biases the interpolation stencils adaptively according to the local smoothness of the data should be employed instead. Hence, to maintain the formal high order of the numerical scheme, a more computational expensive ENO interpolation process (see [15-18, 32, 33, 35, 36] for details) is employed to avoid the interpolation across any potential discontinuity. Furthermore, since the data f (X I ) are distributed uniformly on a two-dimensional Cartesian grid, a one-dimensional ENO interpolation process can be applied to the two-dimensional data in a tensorial, aka, dimension-by-dimension, manner.
Hence, the Lagrange interpolation polynomial can be written as where f km = f (x k , y m ), and l k (x) and l m (y) are the one-dimensional Lagrange interpolation polynomial of degrees K and M, respectively. The index k (m) will be biased according to the local smoothness of the solution that would minimize two (K + 1)th-order ((M + 1)th-order) divided differences of f km if a data point is added either to the right (top) or to the left (bottom) of the K th-order (Mth-order) divided difference stencil, and n i (n j ) represents the left-most (bottom-most) index. The process of moving the indices k and m above will give the adaptive ENO interpolation. This is essentially a one-dimensional process, and hence we will describe it in a generic onedimensional framework by assuming the data is a function of x only and independent of y for simplicity.
We will begin by giving the construction of the K th-order divided difference; • The 0th-order divided difference of f i is defined by: • The kth-order divided difference is then defined recursively by: The interpolation polynomial, in general, can be written as The n i is the left-most index determined by the local smoothness of f (x) at x i as measured by the K th-order divided differences.
• For the starting point n i with the given x 0 : We first locate the nearest neighboring point is formed across the two-points ENO stencil S 2 = {x n i , x n i +1 }.

Numerical Results
To demonstrate the numerical performance of the cp-and tc-embedding methods (refer collectively as embedding methods), we present several classical benchmark problems in one-dimensional scalar hyperbolic conservation laws with a linear, and both convex and non-convex nonlinear fluxes F(Q M ) on a manifold. The tested problems are given on three representative manifolds with an increasing complexity where R x and R y are the radii in the x-and y-directions as given in Table 1.  Table 1 The radius (R x , R y ) with α = 3 2 θ + π 8 , β = e cos(2θ) , the arc-length C k , the computational Cartesian domain x × y, the numbers of grid points N x × N y , and the mesh resolution x = y of the three manifolds M 1 , M 2 , and M 3 The manifolds M 1 , M 2 , and M 3 are a unit circle, an ellipse, and a complex closed curve, respectively, as illustrated in Fig. 2.
The counterclockwise direction is defined as a positive direction. Therefore, the tangential vector at (x, y) is v M (cp(x)) = p (θ ) || p (θ )|| 2 ∈ T M , and the arc-length 1 of the manifold M k is Figure 3 shows the global and local extended velocity field v E ∈ around the manifold M 3 as depicted by (13). The velocity v E increases (decreases) in speed gradually from the v M on the manifold away radially outward (inward) according to the locally signed curvaturẽ κ(cp(x)) and follows the same direction as the direction of the tangent vector at a given point The computational Cartesian domain x × y, the number of grid points N x × N y , and the mesh resolution x = y and other relevant parameters employed in this study are summarized in Table 1 unless specified otherwise.
Under this spatial discretization setting, the corresponding collocating points on the manifolds will not be uniformly distributed and tend to cluster around a large curvature. For presentations of the results in the following numerical experiments, by the second-order ENO interpolation from the collocation points on the manifold, we will obtain and plot 500 data points of the numerical results and the reference solutions distributed evenly along the angle θ of the manifolds. By considering a (local) parameterization γ (s) with s being the arc-length of the manifold The SPDEs of scalar hyperbolic type (1) can be written in the equivalent form of We will begin by examining the accuracy and rate of convergence of the cp-and tcembedding methods in solving a linear wave advection with a smooth solution on the manifolds. We will then examine the essentially non-oscillatory shock-capturing performance by solving the inviscid nonlinear Burgers' equation with a composite-waves initial condition on the manifold M 1 . The numerical dissipation property and the resolution efficiency in resolving fine-scale structures in the complex solution are also studied.
Finally, several sets of benchmark one-dimensional shock-tube problems are simulated on M 3 with the tc-embedding method are presented. In particular, the nonlinear hyperbolic conservation laws with a convex flux (e.g., inviscid Burgers' equation) and a concave flux (e.g., traffic flow problem) that forms isolated shock or isolated rarefaction wave, and a nonconvex flux (e.g., Buckley-Leverett equation) that forms compound wave are used to validate the performance of the tc-embedding method on a complex manifold.

The Accuracy, Resolution, and Efficiency of the Embedding Methods
As for any newly developed numerical scheme, it is essential to examine the accuracy and rate of convergence of the scheme in solving a linear wave equation with a smooth solution. The L 1 , L 2 , and L ∞ errors of the embedding methods in solving the linear wave equation show quantitatively the accuracy and the rate of convergence on the three manifolds M 1 , M 2 , and M 3 . For a non-smooth solution, the nonlinear inviscid Burgers' equation with the composite multi-waves initial condition is solved using the embedding methods to assess their relative performance in essentially non-oscillatory shock-capturing and resolution of smooth finescale structures. Their results are compared qualitatively with the reference solution on the manifold M 1 . The CPU times of the embedding methods are also obtained to demonstrate their relative computational efficiency.

Example 1: Linear Wave Equation
The linear wave equation with a smooth solution is used to obtain the errors in examining the accuracy and the convergence rate of the embedding methods.
The linear flux F is and the initial condition is Q 0 (s) = sin(2πs/C k ).
The exact solution is Q M (s, t) = Q 0 (2π(s − t)/C k ). The final time is T = C k . Figure 4 shows the errors without the nonlinear minmod slope limiter (ω = 0 in (25)) in the CFV scheme since the solution is smooth. In the figure, a thick black line with a slope of two is included for easy reference when comparing the convergence rate of different methods. The errors of the cp-and tc-embedding methods are shown with the dashed and solid lines, respectively. It is easy to see that the errors of the tc-embedding method are about ten times smaller than the errors of the cp-embedding method for all three manifolds. The convergence rates of both embedding methods, as expected, are second-order.
Since we are attempting to solve the hyperbolic conservation laws with non-smooth solutions, the nonlinear minmod slope limiter needs to be activated to suppress any Gibbs oscillations that might appear in the solution. Figure 5 shows the errors in the CFV scheme with the nonlinear minmod slope limiter (ω = 1.25 in (25)) activated. Here, we include another thick black line with a slope of one as the limiter might reduce the convergence rate by a half order. It is easy to see that the errors of the tc-embedding method are still smaller than the errors of the cp-embedding method for all three manifolds. The convergence rate is reduced slightly to around 1.5 in the L 2 norm and 1.2 in the L ∞ norm. These results are  [19]. In summary, the results here show that the tc-embedding method is more accurate than the cp-embedding method, and both methods converge at the same rate as the mesh is refined.

Example 2: Inviscid Burgers' Equation with Composite Waves IC
The inviscid Burgers' equation with the composite waves initial condition is employed to demonstrate the essentially non-oscillatory shock-capturing capability near the discontinuities and the high-resolution at the fine-scale smooth structures of the scheme.
The nonlinear convex flux F is and the initial condition consisting of composite waves is  with η = 4 2s C k − 1 . The final time is T = 1 20 C k . Since there is no exact solution to the Burgers' equation with this complex initial condition, the solution of the cp-embedding method with the mesh resolution x = y = 1 250 is used as the reference solution. In Fig. 6, the solution Q M (s, T ) at the final time T = 1 20 C 1 on the manifold M 1 computed by the embedding methods with the mesh resolution x = y = 1 30 and x = y = 1 40 are shown. The solution on the manifold exhibits multiple shocks and fine-scale structures. Both embedding methods can capture multiple discontinuities in an essentially non-oscillatory manner as the solution evolves in space and time, forming discontinuities from an otherwise smooth solution. If not well-captured by the locally adaptive numerical dissipation at a strong gradient, the Gibbs oscillations will pollute the smooth fine-scale structures. Furthermore, as shown in the right sub-figure of Fig. 6, the tc-embedding method has less numerical dissipation than the cp-embedding method allowing an efficient resolution of the fine-scale structures, which is a very important feature for a long time simulation.
Finally, we present the CPU times of the embedding methods to demonstrate their relative computational efficiency. Both embedding methods are written in Matlab version R2019a and ran on a desktop computer equipped with a four cores Intel(R) Core(TM) i7-6700 CPU with a clock speed of 3.40 GHz. Tables 2 and 3 show the CPU times of the cp-embedding method (CPU cp ) and the tc-embedding method (CPU tc ), and the speedup factors (SF = CPU cp /CPU tc ) for an increasing resolution. The CPU times are obtained by running the Example 1 on the unit circle M 1 to the final time T = C 1 and the complex manifold M 3 with 20 Runge-Kutta time steps. We shall also denote n(X I ) and n(X G ) as the number of cells in the sets X I and X G , respectively. The total CPU times (CPU total ) of the embedding methods are mainly composed of the CPU times used in the CFV method (CPU cfv ) and the ENO interpolation (CPU eno ), that is, CPU total ≈ CPU cfv + CPU eno . In terms of CPU cfv , the embedding methods are basically the same (SF ≈ 1). In terms of CPU eno , the tc-embedding method is about 3.5 times faster than the cp-embedding method (SF ≈ 3.5). Although n(X I )/n(X G ) ≈ 3, the CPU eno dominates over the CPU cfv . In terms of CPU total , the tcembedding method is about 2.5 times faster than the cp-embedding method (SF ≈ 2.5). The reason for such large speedup factors is that the former eliminates the costly computations in reconstructing the CAN-property for all x ∈ X I at each RK-TVD stage required by the latter.
In summary, the tc-embedding method has better accuracy, improved resolution, and reduced CPU time than the cp-embedding method.

The One-Dimensional Benchmark Shocked Flow Problems on M 3
This section presents the numerical results of solving the hyperbolic conservation laws on M 3 with several benchmark one-dimensional scalar hyperbolic conservation laws with discontinuous solutions. Both the linear (F = 0), convex (F > 0), concave (F < 0), and non-convex (F changes sign) fluxes of the hyperbolic conservation laws are considered. The problems are the linear advection equation with a discontinuous initial condition (linear flux), the inviscid Burgers' equation with a stationary shock and with an interaction of the moving shock and composite waves (convex flux), the traffic flow problem (concave flux), and the Buckley-Leverett equation (non-convex flux). These problems demonstrate the robust performance of the tc-embedding method in resolving fine-scale structures efficiently even in the presence of a shock and the essentially non-oscillatory capturing of complex shocks and rarefaction waves existed in the nonlinear PDEs. As the results computed on the other two manifolds M 1 and M 2 are essentially the same as the one computed on the M 3 when the final time T is scaled proportionally with their respective arc-length C k ; they are omitted here for saving space. In

Example 3: Linear Advection Equation with Multi-Waves IC
In this example, the simple linear advection problem is solved with the linear flux F(Q M ) and a piece-wise continuous multi-waves initial condition Q 0 (s) in (33) as given below.
The linear flux F is and the multi-wave initial condition consists of a smooth Gaussian function, a discontinuous square function, a piece-wise continuous triangle function, and a smooth elliptic function, with a = 0.5, z = −0.7, δ = 0.005, α = 10 and β = log 2 36δ 2 . The exact solution is Q M (s, t) = Q 0 (2(s − t)/C 3 − 1). The final time is T = C 3 .
As shown in Fig. 7, the tc-embedding method is capable of advecting the discontinuous solution in an essentially non-oscillatory manner while resolving the piece-wise smooth solutions well spatially and temporally. The result is very encouraging.
Next, we shall solve the SPDEs with a nonlinear convex and non-convex flux, giving rise to complex finite-time singularity even with a smooth initial condition.

Example 4: Inviscid Burgers' Equation with Stationary Shock
Inviscid Burgers' equation is the simplest prototype for the scalar hyperbolic conservation laws with a nonlinear convex flux that the solution can develop either an isolated shock wave (discontinuity in the solution) or an isolated rarefaction wave (discontinuity in the first derivative of the solution). They are the weak solutions of the inviscid Burgers' equation and satisfy an entropy inequality with entropy function η = Q 2 [20]. The Burgers' equation is a fundamental partial differential equation in applied mathematics, such as fluid mechanics and gas dynamics. Here we solve (33) on M 3 with the nonlinear convex flux F and a smooth initial condition Q 0 (s) given below with the tc-embedding method.
The nonlinear convex flux F is and the smooth initial condition is The final time is T = 1 4 C 3 . Figure 8 shows the numerical results of the simulations. As shown in the space-time evolution of the solution Q M (s, t) in the middle sub-figure, the initial smooth function develops into a stationary discontinuity (shock) at arc-length s = C 3 /2 and the wave breaking time t = C 3 /(2π), as predicted by the theory when the left-and right-going characteristics have collided together. The weak solution also satisfies the Rankine-Hugoniot condition for the shock speed. The simple piece-wise continuous numerical solution has no discernible Gibbs oscillations. The example demonstrates that the tc-embedding method can capture the shock wave adaptively in an essentially non-oscillatory manner while resolving the piecewise smooth solutions well spatially and temporally.
This problem is relatively easy to handle as the discontinuity is stationary in space and time. Next, we shall tackle the same problem but with moving shocks and with a complex initial condition consists of composite waves. This problem will challenge the space-time adaptivity characteristic of the tc-embedding method.

Example 5: Inviscid Burgers' Equation with Composite Waves IC
This example is given to show the performance of the tc-embedding method on M 3 . The initial condition is complex, including four discontinuous with different jump ratios and two smooth sine functions with different wave frequencies.
The nonlinear convex flux F is and the initial condition consisting of composite waves is given in (37). The final time is T = 1 20 C 3 . Figure 9 shows the Burgers' equation with the composite waves initial condition. As shown in the space-time evolution of the solution Q M (s, t) in the middle sub-figure, the left-most sine wave develops into stationary shock at s = C 3 /8, while the right-most high-frequency sine wave to the right of the strong shock develops into a train of shocklets. The two rightmoving shocks are formed in the middle and moves to the right. They are captured essentially non-oscillatory as they evolve around the manifold M 3 in space and time, as shown in the right sub-figure. This challenging example demonstrates the space-time adaptivity characteristic of the tc-embedding method on a non-trivial shape of the manifold.

Example 6: Traffic Flow Problem
In mathematics and transportation engineering, traffic flow is a fundamental theory to analyze and study the interactions between travellers (mainly cars) and infrastructure, with the aim of understanding the relationship between traffic flow, velocity and density, reducing the occurrence of accidents, and improving the efficiency movement of traffic.
The nonlinear concave flux F is where Q M is the cars' density. We take v max = 1 and Q max = 1 as the limited speed and density, respectively. The initial condition is where η = 2s C 3 − 1. The final time is T = 2 5 C 3 . As shown in Fig. 10, there are two strong discontinuities (front and back of the traffic jam) in the initial conditions. The left discontinuity (shock) propagates to the left as the cars in the low-density region slow down and become congested as cars enter the high-density region. The right discontinuity evolves into a rarefaction wave as cars' density reduces gradually, as the cars speeding up, from the high-density region toward the low-density region. The left-going shock and the left-going rarefaction waves are well-captured in an essentially nonoscillatory manner by the tc-embedding method on M 3 spatially and temporally, as shown in the middle and right sub-figures.
The examples above show that the tc-embedding method performs equally well for both the convex and concave fluxes. We will conclude the study showing that the tc-embedding method also works well with a non-convex flux, where the second derivative of the flux changes sign, on M 3 .

Example 7: Buckley-Leverett Problem
In fluid dynamics, the Buckley-Leverett equation is a model for a two-phase oil-water immiscible displacement process in porous media in a one-dimensional or quasi-one-dimensional reservoir.
The nonlinear non-convex flux F is with the viscosity μ = 1/4. The initial condition is where η = 2s C 3 − 1. The final time is T = 1 5 C 3 . The flux is S-shaped with an inflection point and, therefore, non-convex. It yields a rarefaction-shock compound wave, also known as the Buckley-Leverett (BL) profile, which consists of a shock wave immediately followed by a rarefaction wave (see [20] for details). Figure 11 shows that the BL profile is captured essentially non-oscillatory and sharply by the tc-embedding method. The space-time evolution of the solution shows the formation process of the BL profile as the shock waves move towards the right stably together with the expanding rarefaction fans attached.

The Burgers' Equation on the Two-Dimensional Manifolds
This section presents preliminary numerical results for applying the tc-embedding method for solving the moving Burgers' equation on two-dimensional manifolds embedded in R 3 . More specifically, they are the torus-shaped and spherical-shaped manifolds. We conjuncture that the generalization of Theorem 2 to R 3 can be done by using extended velocity field as in (13) on codimension-2 manifolds constructed by slicing M along v M direction. We leave the proof to our future work.
We first illustrate the performance of the tc-embedding method in solving SPDEs on a two-dimensional parametrized torus surface: x(θ, φ) = (R + r cos φ) cos θ, y(θ, φ) = (R + r cos φ) sin θ, z(θ, φ) = r sin θ. (47) where the parametrization parameters are θ ∈ [−π, π) and φ ∈ [−π, π). R = 1 is the distance from the center of the tube to the center of the torus, and r = 0.5 is the radius of the tube. The SPDEs Burgers' equation simulates the finite-time formation of shock waves and their movement on a surface. The initial condition is The final time is T = π. In Fig. 12, the solutions Q M (x, t) of the Burgers' equation at times (t = 0, 1, 2) on the torus surface computed by the tc-embedding methods with the mesh resolution x = y = z = 1 25 are shown. The results show that the tc-embedding method can capture the formation of the strong gradient and shock sharply and essentially non-oscillatory and its later evolution on the torus-shaped manifold.
(49) where the parametrization parameters are θ ∈ [0, 2π) and φ ∈ [0, π). R = 1 is the radius. The initial condition is where C = 10 −5 , a 1 = π 4 , k i = 2i − 1, a 2 = π 5 , h i = i, and σ = π 8 . The final time is T = π. In Fig. 13, the solutions Q M (x, t) of the Burgers' equation at times (t = 0, 1, 2) on the spherical surface computed by the tc-embedding methods with the mesh resolution x = y = z = 1 25 are shown. The initial condition (t = 0) shows the four peaks on the spherical surface. Similar to the single shock formation in the previous example, the proposed algorithm can capture the formation of the strong gradient and shock sharply and essentially non-oscillatory. The movements and mutual interactions among the Gaussian functions are well-resolved, and discontinuities are well-captured, as depicted in the middle and right subfigures.

Conclusion
A new closest point time-continuous embedding (tc-embedding) method is developed for solving nonlinear scalar hyperbolic conservation law (PDEs) on one-dimensional, connected, smooth, and closed manifolds (SPDEs). It improves upon the classical closest point embedding (cp-embedding) method that requires a costly re-establishment of the exten-sion function's constant-along-normal (CAN) property at every time step; the tc-embedding method incorporates the CAN-property analytically and explicitly in the SPDEs directly. The main idea in the tc-embedding method for solving SPDEs on a one-dimensional manifold is that under certain minor restrictions, the time-independent extended velocity field only needs to be determined once at the initialization stage and then used subsequently for all times in the method; thus, a substantial saving in the computational times can be realized.
Due to the nonlinear nature of the PDEs, finite-time singularities, such as shocks and rarefaction waves, can appear spontaneously in the solution. A second-order central finite volume scheme with a nonlinear minmod slope limiter is used to approximate the flux gradients discretely in a computational tube placed inside a Cartesian rectangular uniformly spaced mesh to mitigate Gibbs oscillations. It requires values at the ghost cells to compute the flux gradients at the boundary cell at the edges of the computational tube. By using the cp-mapping and CAN-property to map the location of a ghost cell to a point on the manifold, the ENO interpolation process updates the values of the ghost cells using data around the point on the manifold. The adaptive ENO interpolation is an essential technique in obtaining an essentially non-oscillatory value by avoiding interpolating across any potential discontinuity or high gradient in the solution. The third-order total variation diminishing Runge-Kutta (RK-TVD) scheme integrates the numerical solution in time.
We compared both methods' accuracy on the linear advection equation with a smooth solution on simple and complex manifolds. While both methods are formally second order, the tc-embedding method is more accurate than the cp-embedding method overall. We have also compared both embedding methods' resolution power and CPU times in solving the Burgers' equation with a piece-wise continuous solution on a circular manifold. The results indicate that the tc-embedding method has a lesser numerical dissipation, resolves the finescale structures more efficiently, and is roughly two and half times faster CPU times than the cp-embedding method. Finally, several one-dimensional benchmark shocked flow problems, such as the Burgers' equation with convex flux, the traffic flow problem with concave flux, and the Buckley-Leverett equation with non-convex flux, are used to demonstrate the essentially non-oscillatory shock-capturing property and the efficient resolution of fine-scale structures of the tc-embedding method on a complex manifold. The numerical results of the Burgers' equation are also given for the two-dimensional torus-shaped and spherical-shaped manifolds.
These preliminary results show the potential synergy between the embedding method for solving SPDEs with smooth solutions and the nonlinear shock-capturing scheme for hyperbolic conservation laws with non-smooth (discontinuous) solutions. This study serves as the starting point, and there are more challenges in the theoretical and numerical aspects of combining both for solving more sophisticated hyperbolic SPDEs. As immediate follow-up work, extending the theoretical and numerical frameworks of the tc-embedding method for solving the scalar hyperbolic conservation laws on a two-dimensional complex connected manifold with a well-defined velocity field is challenging and will be investigated. As for future works, we will study if the tc-embedding method can help solve the one-dimensional or multi-dimensional system of hyperbolic conservation laws, for example, shallow water equations, Euler equations, multi-phases flow, and multi-components flow, on a stationary or moving manifold. Moreover, high-order shock-capturing schemes, with improved resolution, for example, the weighted ENO (WENO) finite volume/finite difference schemes [5,6,34] on a regular Cartesian mesh and the discontinuous Galerkin [11] method (DG) on a structured and unstructured triangular mesh, is also currently under consideration. support of this research by the National Natural Science Foundation of China (11871443) and a Hong Kong Research Grant Council GRF Grant. The author (Don) also likes to thank the Ocean University of China for providing the startup funding (201712011) to support this work. The authors also thank Mr. Kang-Bo Tian for sharing the program on the tc-embedding methods on the two-dimensional manifold in R 3 .

Author Contributions
The authors have contributed equally in this scholarly work.

Funding
The authors have not disclosed any funding.

Data Availability
No data sets were generated or analyzed during the current study.

Conflict of interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.