A Well-Balanced SPH-ALE Scheme for Shallow Water Applications

In this work, a new discretization of the source term of the shallow water equations with non-flat bottom geometry is proposed to obtain a well-balanced scheme. A Smoothed Particle Hydrodynamics Arbitrary Lagrangian-Eulerian formulation based on Riemann solvers is presented to solve the SWE. Moving-Least Squares approximations are used to compute high-order reconstructions of the numerical fluxes and, stability is achieved using the a posteriori MOOD paradigm. Several benchmark 1D and 2D numerical problems are considered to test and validate the properties and behavior of the presented schemes.


Introduction
The shallow water equations (SWE) with non-flat bottom geometry problems suppose a hyperbolic system of conservation laws with a source term, also called balance laws. In these problems we can achieve steady state solutions in which the flux gradients are not zero, but rather they are exactly cancelled by the contribution of the source term. In practice, the exact cancellation of the flux gradients with the source term is difficult to achieve numerically, if there are no additional considerations when it comes to discretize the source term [1][2][3][4][5][6][7].
When a non-well-balanced scheme is used in SWE with variable bathymetry some fixedvalued artificial perturbations will appear throughout the domain [8]. These perturbations are proportional to the mesh size or to the distance between particles in meshless methods. Well-balanced schemes [9] or schemes that satisfy the C-property, as they were called in early works [10,11], are those which can successfully compensate (balance) the source term with the flux gradients and reproduce steady state solutions. Hence, a well-balanced scheme should preserve static states over time regardless of the bathymetry. Balancing the SWE is a

Governing Equations
The SWE with non-flat bottom geometry suppose a hyperbolic system of balance laws defined as follows: where U U U is the vector of conservative variables and F F F = (F F F x (U U U ), F F F y (U U U )) is the flux tensor, that are defined as The source term vector, S S S(U U U ), accounts for the effects of the bottom bathymetry and is defined as The terms h and g stand for the water height and the gravity acceleration. The position is defined by r r r = (x, y) T , and the elevation of the bottom bathymetry is defined by a smooth function b(r r r ) that is included in the source term. The water flow velocity vector is u u u = (u, v) T , and the elevation of water surface is denoted as η = h + b. In Figure 1 a schematic diagram with h(r r r ), b(r r r ) and η(r r r ) is shown.

A Well-Balanced, High-Order SPH-ALE Scheme
In SPH-ALE methods the problem's domain is discretized into a set of particles that hold the values of the variables in their positions. We can solve the SWE described in section 2 for a given particle i as the weighted sum of the interactions between this particle i and the neighbouring particles j within its influence domain, denoted as stencil. This interaction between particles is conceived as a Riemann problem (RP) at the integration point i j at midpoint between particles i and j. In Figure 2 a schematic representation of the fluid particles, the particle i, its stencil, the neighbour particle j, and the integration point i j located at midpoint between i and j are shown. The positions of these points are denoted by r r r i , r r r j and r r r i j respectively. The normal vector of the i j path is noted as n n n i j = r r r j −r r r i ||r r r j −r r r i || . In the present work, we use the SPH-ALE formulation [17][18][19] with the correction term proposed in [34]: dr r r i dt = where r r r i is the position vector of particle i, G G G i j is the numerical flux tensor computed by a Riemann solver in the integration point i j, which will be defined later, and is the Lagrangian flux tensor of the SWE in a reference moving frame with velocity w w w. It is defined as The terms w w w i , w w w j refer to the Lagrangian velocity of particles i, j, while the Lagrangian velocity at the integration point i j is w w w i j , which will be defined in section 3.1. In an Eulerian frame w w w i = w w w j = 0, and in a Lagrangian frame they are coincident with particle velocities u u u i ,u u u j . The flux tensor defined in Eq. (5) leads to a flux difference formulation of the SPH-ALE scheme [34]. In Eqs. (2), (3) and (4), the term V j represent the volume associated to the particle j, and W i j is a kernel function. Here, the kernel function is defined as In Eq. (6), D is the number of space dimensions. The value of the kernel depends on the positions of the particle i and its neighbouring particles j, as q i j = ||r r r j − r r r i || l i j . The value of l i j is defined from the smoothing lengths associated to the particle i and particle j and related to their volume as The value of the constant C depends on the number of spatial dimensions. We use C = 1 for 1D problems and C = 10 7π for 2D problems. We refer to [24] and its references for a more exhaustive and complete description of the SPH kernel formulations used in this work.

Discretization of the Numerical Fluxes
The numerical flux tensor G G G i j is computed by a Riemann solver, considering that the left and right states of the RP at the integration point i j are given by the polynomial reconstructions of the variables at the integration point i j. In the following, for a given variable α, we consider α − i j and α + i j as the left and right states of the RP respectively. In Fig. 3 a schematic diagram of the mentioned RP and the reconstructions at i j is shown. In this work, the Rusanov flux is considered [35].
In Eq. (7) the term S * i j is the maximum eigenvalue of the flux's Jacobian matrix [35], which in the ALE framework reads as where w w w i j = 1 2 w w w + i j + w w w − i j is the Lagrangian velocity associated to the integration point, c ± i j = g · h ± i j denotes the left and right states of the wave celerities in the RP. Moreover, are the flux approximations of H H H on the left and right sides of the integration point (located at the midpoint between particles i and j), r r r i j with the positive orientation given by r r r j − r r r i and represented by the normalized vector n n n i j , while U U U + i j − U U U − i j is the jump of the reconstructed conservative variables.

High Order Reconstructions with Moving Least Squares
The left and right states at integration points of the conservative variables (U U U − i j and U U U + i j ) and the bathymetry (b − i j and b + i j ) are computed using a cubic Taylor polynomial. These polynomial reconstructions are used to calculate the fluxes H H H − i j and H H H + i j in equation (7), as well as the source term (that will be described in the next section).
The derivatives required to compute the Taylor polynomials are obtained using MLS approximations. In the following, we present very succinctly the MLS approximations approach, and we refer the reader to [21,[23][24][25][26] for a complete description of the method.
We can approximate the value of the variables U U U at a given point, r r r = (x, y) T , using MLS reconstruction as where N i is the number of neighbour particles in the stencil of particle i and φ j (r r r ) are each of the MLS shape functions associated o the neighbour particles j [21,25]. These shape functions, are calculated in function of the relative positions of the particles j and i and they are weighted by a kernel function. The shape functions, gathered in vector where W W W D (r r r ) is a diagonal matrix which is obtained from the kernel function evaluated at r r r j − r r r i for the N i neighbouring particles. In the numerical applications of the present work, we use the truncated exponential kernel function [23,24] W i j = W (r r r i , r r r j ) = e 1−β 2 − 1 e − 1 ; β = |r r r j − r r r i | 2 · max(|r r r k − r r r i |) , ∀k ∈ N i

Well Balanced Methods: Source Term Discretization
In order to obtain a well balanced scheme, a new discretization of the source term is presented. The proposed discretization is based on the splitting of the source term, and it adapts to the SPH-ALE scheme the source term separation procedure presented by Xing and Shu for FD methods [7] and to FV and DG methods [6]. If we consider the contribution of the non-flat bathymetry and no other external forces, such as bed friction, the source term on a 2D SWE scheme reads: The discretization of S 2 and S 3 components of the vector (10) is the same, so in the following only the S 2 term is derived. First, we proceed to separate the source term into two different terms: Usual kernel approximations of a function and its gradient, read as In order to ensure zeroth-order consistency, it is usual to approximate the gradient as [17,36].
The value of a function at the midpoint between particles can be computed as and then, Thus, applying this approximation to the gradient of b, ∇b = ( ∂b ∂ x , ∂b ∂ y ) T , the term S A 2 , only involving the first component of the gradient, can be expressed as: An analogous procedure applied to the S B 2 term leads to: If we add both split terms (12,13) into a single sum, we finally obtain S 2 , and S 3 as the proposed discretization of the source term: Applying the following discretization of the source term to the SWE hyperbolic system (7), the flux gradient is explicitly canceled by the source term when water at rest conditions are considered (h + b = η =constant; u = v = 0 for all the particles). This is strictly true for the equations derived from the conservation of momentum (those that consider S 2 and S 3 respectively). The equation of the water height, derived from mass conservation, is examined in the following. We consider the water height equation, and we define and the high-order reconstructions Then, the water height equation reads as Considering water at rest conditions, we obtain: Note that w w w i j = 0 0 0 and w w w i = 0 0 0 for the water at rest conditions for both Lagrangian and Eulerian approaches. Eq. (16) shows that the difference of the jump term is not exactly neglected for the water at rest conditions, as the reconstructed values of the water height (h) may be different in the presence of non-flat bottom bathymetry, and the resulting numerical scheme is not well-balanced. In this work, following the works presented in [5,37], we propose a well-balanced SPH-ALE discretization of the mass conservation equation. We consider a non-erodible bathymetry hypothesis, and under this hypothesis, the transport operators of the water height (h) and the bathymetry (b) can be expressed as Then, the mass conservation equation in terms of h reads and, introducing the operator L w (h) and the flux Since η = h + b and the bathymetry is non-erodible we can rewrite the mass conservation equation in terms of η as We can discretize Eq. (17) in particle i approximating the ALE flux term as a numerical flux [17] including the kernel gradient correction [34] as follows Note that we use the Rusanov numerical flux [35] considering the water surface elevation η as the updating variable on time, since Eq. (17) is written in terms of this variable.
We can now revert the change of variables to write Eq. (18) in terms of water height h before proceeding with the time discretization Using that Finally, substituting in Eq. (18) we get When water at rest conditions are considered, we obtain Hence, since η is constant, an exact zero mass flux is granted. Thus, we have derived a new SPH-ALE method in which the numerical flux for the mass conservation equation is discretized considering the water surface elevation η as the updated variable, instead of h. It is important to note that in the proposed method the conservative variable is h. This is different of the procedure followed by the pre-balanced methods [38][39][40][41] in which the conservative variable is changed for the water elevation η.

Final Discretization
To summarise, and for the sake of clarity, we write down the system of equations to be solved and discretizations adopted in the presented method.
And the source term is discretized as

Time Discretization
A TVD third-order Runge-Kutta time discretization [42] is used in this work where F (α) is the spatial operator, and α is any of the time dependant variables from Eq. (19). In this work the time step is defined according to

MOOD Method: A posteriori stabilization
High order reconstruction methods, such as the fourth-order MLS polynomial reconstruction employed in this method, may produce artificial oscillations in the vicinity of shocks or discontinuities of the solution. In order to avoid these non-physical oscillations, some kind of stabilization procedures must be introduced. The classic approaches for stabilization are based on artificial viscosity terms or slope limiters, that are estimated a priori. This means that the method predicts at time step n the cells/particles/elements that will become problematic in the next time step n + 1.
In the present work, the a posteriori paradigm is used. In particular, we use the multidimensional Optimal Order Detection (MOOD) limiting procedure [24,32,33,[43][44][45][46]. In this a posteriori method, a candidate solution is obtained by solving the system of equations using the most accurate method at disposal without any stabilization mechanism.
This candidate solution may contain N a N values or excessive numerical oscillations. In order to prevent that, some detection chains are applied to that candidate solution to identify the troubled particles.
After the MOOD detecting procedure, the fluxes are recalculated considering a lower reconstruction order for those particles which were detected. This procedure (MOOD loop) continues until all the particles of the candidate solution successfully pass the detection chain. It has been shown that the schemes built with this approach are convergent and stable, provided that the last scheme used is convergent and stable. Thus, the last scheme (termed as the parachute scheme [32]), must be convergent and stable. In our case, the parachute or base scheme is the first-order SPH-ALE method. In the present implementation of the MOOD procedure, if a particle is marked as problematic, the order of the reconstructions switches from fourth to first order.
In this work, we implemented a MOOD loop based on the method proposed in [24]. The first detector of the detection chain is the Physical Admissible Detector (PAD). The PAD ensures the positivity of the water height and of the particle volume, and also search for N a N values. The second detector is the Discrete Maximum Principle (DMP) that detects local maximum of the water elevation (η) as a potential numerical oscillation. Finally, the Plateau Detector (PD) relax the DMP detector by admitting oscillations under a margin of 10 −15 , that is, machine truncating errors that would mistakenly trigger the DMP detector. For a comprehensive explanation of the implementation of the MOOD method in the SPH-ALE scheme, we refer the interested reader to [24].

Bathymetry Reconstruction
When a Lagrangian approach is considered, the location of the particles changes over time and thus, the points where the bottom bathymetry must be evaluated are different each time step. At the same time, in practical applications of the SWE, the bottom bathymetry is only known at a discrete set of points. Hence, a SWE method must consider a procedure to obtain the bathymetry at any given point inside the problem's domain.
In order to design an algorithm able to deal with scattered bathymetry data, we propose a new procedure to reconstruct the bathymetry from a discrete set of points. We consider that the bathymetry is only known at a discrete set of points that we call bathymetry points. In our numerical tests the bathymetry points are coincident with the initial location of the fluid particles but, any other locations are also valid for the proposed algorithm. We approximate the bathymetry at any given position r r r i inside of the problem's domain employing m th -order Taylor polynomial reconstructions, namely b p (r r r ) where the gradients and derivatives are computed with MLS from its closest bathymetry point r r r p b(r r r where ∇ n b 0 p (r r r ) are the MLS derivatives of the bathymetry at r r r p , b 0 k are the values of the bathymetry at the bathymetry points r r r k , and Φ 0 k are the MLS shape functions of each neighbouring bathymetry point of r r r p , M p . The m th -order polynomial reconstruction of the bathymetry based on r r r p is denoted as b p (r r r ). Note that we use this procedure even though in the proposed numerical examples, an analytical expression for the bathymetry is known.
Then, a MOOD procedure for the bathymetry approximation starts with a candidate bathymetry computed using a fourth-order polynomial reconstruction. Then, the DMP and PD detectors are applied to assess that the approximation is valid in terms of smoothness. If the approximation is identified as troubled by the detectors, a new candidate approximation employing a lower order MLS reconstruction is checked. This MOOD loop continues until the candidate solution successfully passes the detection chain. This MOOD procedure grants a smooth bathymetry reconstruction, compatible with the formulation presented in Sect. 2.

Numerical Results
This section presents the numerical results for several benchmark problems which aim to assess the accuracy and efficiency of the proposed method. These problems also aim to assess the Well Balancedness of the method for both, the high-order reconstruction scheme and also the first-order base scheme. Thus, in the following test cases we are employing and comparing the performance of the proposed method, labelled as SPH-ALE-MOOD, with the first-order SPH-ALE scheme which is labelled as the Base scheme.
A Lagrangian approach and homogeneous boundary conditions are used for the proposed test cases, unless indicated in the text. The boundary conditions in this work were implemented using the mirroring ghost particles procedure, as in [47,48].

Convergence Test: Steady Subcritical Flow Over a Bump
In this test the order of convergence of our scheme is checked for a smooth solution. This test case presents a steady subcritical flow over a bump with a smooth bathymetry and it is commonly use to analyze the accuracy and convergence of the numerical methods [30,49,53]. The bathymetry is given by Considering a steady state, we can obtain the analytical solution of the problem as This problem was computed until convergence with different refinements, using a SPH-ALE scheme with fourth-order MLS reconstructions.
In table 1, the L 2 − norm of h and u and the obtained order of convergence are presented for several particle discretizations. In the first set of particles, until 160, tends to fourth-order convergence, but after it decays to the expected second-order [17]. This behaviour is due the fourth-order reconstruction of the variables and the second-order spatial integration of the scheme.

Lake at Rest Test Cases
This set of cases is widely used to test the well balancedness of numerical methods [1,29,50]. These problems test the ability of the schemes to preserve the water at rest solution (null velocity and constant water surface elevation) over time for an arbitrary bottom bathymetry .
In the following test cases, the proposed method is compared with an SPH-ALE method with a non-well balanced discretization of its source term. This non-WB source term is discretized as

1D Lake at Rest Test Cases
In this first test case, we tested two different bottom bathymetries. The first bathymetry b 1 (x), consists of a single bump [1,29,51,52]. The second bottom bathymetry b 2 (x), is more complex and considers a variable ground elevation [53].
The two bottom bathymetries are defined as follows: For these cases, 200 particles were equally distributed throughout the computational domain. All of the particles were given the same water surface elevation η = 1 and null velocity as initial conditions. The CFL number was set to C F L = 0.9. The results of the two bathymetries test cases for T = 1 are shown in Figs. 4 and 5 . As it is shown in the mentioned figures, the non-well-balanced methods introduce artificial perturbations. The high order non-well-balanced method reduce the magnitude of the oscillations, but it still fails to obtain an exact flat surface. The proposed well-balanced SPH-ALE method is able to preserve the stationary static solution over time for both, the high-order and the base schemes, keeping the relative error of water surface elevation at the order of 10 −16 , which is at the level of the machine precision.

2D Lake at Rest Over a Bump
We consider a test case that takes the lake at rest conditions for a 2D domain. This test case involves a circular domain with a radius of R = 45 discretized with 60429 scattered particles. All particles were initialized with lake at rest conditions with η = 1. The CFL number was set to C F L = 0.9. The bottom bathymetry is defined as follows: We run the simulation until a final time of T = 1. Figure 6 shows that the proposed method is able to preserve the stationary static solution over time for both, the base scheme and the SPH-ALE-MOOD methods. Note that all the particles are represented in the figure. The relative error of water surface elevation is of the order of 10 −16 .

Propagation of Small Water Surface Perturbations on Water at Rest Cases
The goal of these problems is to asses the Well Balancedness of the presented method when a small water perturbation is present on the water at rest initial conditions. We consider initial null velocity and constant water surface elevation for all the particles, except in a subregion of the domain, in which a small perturbation is added to the water surface as initial conditions.
The setup of the following test cases were taken from [1]. These test cases have been widely used by other researchers to test different mesh-based numerical methods, in 1D [6,49,54] and in 2D [6,12,54,55]. For SPH methods, these test cases have been used in [29,31]. A C F L = 0.15 was employed in these cases.

Propagation of Small Water Surface Perturbations on 1D Lake at Rest Cases
In this test case, we consider a domain Ω : Throughout the problem's domain, we set 200 particles with null velocity and the following water height  Fig. 7. In Fig. 7 our results are compared with a reference solution obtained by LeVeque [1] employing a finite volume method with a very fine grid. For the small perturbation case, we also plot the results obtained in [29] using a well-balanced SPH scheme with 400 particles. Neither of the proposed methods show any artificial perturbation due to the bathymetry bump, in consonance with the well-balanced behaviour of the schemes. The SPH-ALE-MOOD method obtains a very accurate representation of the solution. It is also important to remark that the MOOD method is employing two different reconstruction orders in different particles simultaneously (fourth-order reconstruction in non-troubled cells and first-order reconstruction in troubled cells) throughout the domain. Despite this, the method is able to keep the well-balanced behaviour. For the small perturbation case, we note that the result   [29] of the proposed SPH-ALE-MOOD method is comparable to that of [29], using half of the particles.

Propagation of Small Water Surface Perturbations on 2D Lake at Rest Case
In this test case, we consider a rectangular domain Ω : (x, y) ∈ [0, 2]×[0, 1] and the bottom bathymetry reads as The computational domain is discretized with an equally distributed set of 200×100 particles, with initial null velocity and the water height defined as The numerical results for T = 0.6, T = 1.2 and T = 1.8 are presented in Fig. 8 both for the SPH-ALE Base scheme and for the SPH-ALE-MOOD methods. In Fig. 9 a 3D view of the results from the SPH-ALE-MOOD method are shown. Figures 8 and 9 show the results obtained using the base scheme and the SPH-ALE-MOOD methods. We notice that both methods obtains a solution free of spurious oscillations, showing the well-balanced character of the proposed schemes. Moreover, the proposed SPH-ALE-MOOD method obtains a very accurate solution, which is comparable to those obtained in the literature [29,55].

Moving Fluid with Height Discontinuity over Constant Bed Slope
This academic numerical test was set to check the performance of the proposed method in an extreme situation when water height discontinuities, fluid velocity and, non-flat bathymetry are present at the same time. Here, we examine if the scheme preserves the well-balance property when shocks, contact discontinuities and, rarefactions are present in the solution. This test case was first proposed in [57] as a modified version of a test from [56]. The initial conditions are defined as The bottom bathymetry is given by The discretization is performed using 240 particles, initially distributed in an extended domain Ω 0 : x ∈ [0, 1.2]. We did this in order to get 200 particles on the computational domain (Ω : x ∈ [0, 1]) at the final computing time T = 0.1, as it is shown in Fig. 10.The CFL number was set to C F L = 0.6. In Fig. 10 the results of this problem are shown for the SPH-ALE-MOOD method compared to reference solution from [57]. It is observed that the proposed method captures accurately the position of the shock. The MOOD procedure is able to prevent the appearance of non-physical oscillations at the surroundings of the shock.

2D Circular Dam Break on Flat Bathymetry
Dam break tests are a classic way to test the accuracy of SWE methods and specially, their ability to correctly solve problems with discontinuities and shock waves. This test case was introduced in [58] and has been used in many different works since then, i.e. [12,59,60]. This problem considers a squared domain [−25, 25] × [ −25, 25] in which we distribute 40090 particles with an scattered disposition. All of the particles were initialized with null velocity, null bottom bathymetry, and the following water height as initial condition: The final results for different times are shown in Figs. 11 and 12 . The CFL number was set to C F L = 0.90. In Fig. 11 the results obtained with the SPH-ALE Base scheme and SPH-ALE-MOOD methods for T = 0.69 are compared with [100×100] FV-HLLC first-order and second-order-MUSCL reference solutions from [60], respectively. We can see that the SPH-ALE-MOOD method's results agree well with the HLLC-MUSCL FV method reference solution. In Fig. 12, the 3D view of the water surface elevation obtained with the SPH-ALE-MOOD method for six different times T = 0.0, 0.2, 0.4, 0.6, 0.8, 1.0 is shown.

2D Circular Dam Break Problem on a Non-Flat Bed
This 2D dam break benchmark problem was firstly introduced in [61], and it was later reproduced in different SWE works such as [12,59,62]. This test case considers a [0, 2]×[0, 2] domain in which 201×201 particles are equally distributed along both axis. All of the particles were given null velocity and the following water surface elevation as initial condition: And the bathymetry is given by: The results obtained with the SPH-ALE-MOOD method are compared in Fig. 13 with a reference solution (fourth-order FV 800 × 800 grid method) digitized from [59]. We see that the presented method provides accurate solutions, specially at detecting the position and strength of the shock wave.

Steady Flow Over a Bump Test Cases
The following cases were first proposed by Goutal and Maurel in [51] and, since then, they have become classic benchmark test problems to test the behaviour of new numerical schemes for the SWE [3,6,30,63].
For this test case we consider a purely Eulerian approach, using 200 particles which are equally distributed in the domain Ω : x ∈ [0, 25]. The gravity acceleration is g = 9.812. The CFL number was set to C F L = 0.9. The bathymetry of the problems is defined as follows  We consider three different flow cases using the same bottom bathymetry: A subcritical flow (SF), a Transcritical flow without a shock (TF), and a Transcritical flow with a shock(TFsh). The Dirichlet boundary conditions and initial conditions for each of the three different flow cases are shown in Table 2. Note that the boundary conditions that are not determined in the table for the inlet and the outlet are considered as null Neumann boundary conditions. The subcritical flow case (SF) remains in the subcritical regime all along the problem's domain. However, in the steady solution of the transcritical flow without a shock case (TF) the flow regime changes from subcritical at the inlet, to supercritical downstream the bump. In the steady solution of the transcritical flow with a shock (TFsh), the flow regime changes from subcritical at the inlet, to supercritical downstream the bump and then abruptly changes again into subcritical thorugh a shock wave. The results of the calculations are presented in Fig. 15. As we can see, the presented SPH-ALE-MOOD method obtains solutions in very close agreement with the analytical solution. The method is able to reach and preserve the steady state solution over time. No artificial perturbations due to the bottom bathymetry are observed, in consonance with the well balanced behaviour of the proposed scheme. The presented solutions are comparable to those obtained using mesh-based methods [2,3,6,53,55] or SPH methods like [30].

Conclusions
A new well balanced SPH-ALE-MOOD scheme to solve the SWE is presented in this work. The presented numerical method is accurate and provides well-balanced solutions. The wellbalanced behavior is obtained through a new discretization of the source term. The proposed discretization of the source term is able to exactly cancel the contribution of the flux gradient for water at rest conditions.
A series of 1D and 2D benchmark problems were considered to assess the performance of the proposed method. It is shown that the proposed methodology is able to solve the SWE with variable bottom bathymetry and discontinuities in the solution. Moreover, the ALE formulation allows to consider Lagrangian and Eulerian approaches to solve the SWE, and provides high-accurate, well-balanced and robust solutions for both purely Eulerian and Lagrangian approaches. The MOOD limiting procedure is able to prevent artificial oscillations in the vicinity of discontinuities and shocks, preserving the well-balanced behaviour of the proposed scheme.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. This work has been partially supported by FEDER funds of the European Union, Grant #RTI2018-093366-B-I00 of the Ministerio de Ciencia, Innovación y Universidades of the Spanish Government and by the Consellería de Educación e Ordenación Universitaria of the Xunta de Galicia (grant# ED431C 2018/41). Luis Ramírez also acknowledges the funding provided by the Xunta de Galicia through the program Axudas para a mellora, creación, recoñecemento e estruturación de agrupacións estratéxicas do Sistema universitario de Galicia (reference # ED431E 2018/11).

Availability of Data and Material
The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.

Conflict of interest
The authors declare that they have no conflict of interest.
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. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.