* Peñaloza et al. (2018) Peñaloza, C. H., Clark, P. C., Glover, S. C. O., and Klessen, R. S. (2018).  CO line ratios in molecular clouds: the impact of environment.  _MNRAS_ 475, 1508–1520.  10.1093/mnras/stx3263
  * Pelupessy and Papadopoulos (2009) Pelupessy, F. I. and Papadopoulos, P. P. (2009).  Molecular Gas, CO, and Star Formation in Galaxies: Emergent Empirical Relations, Feedback, and the Evolution of Very Gas-Rich Systems.  _ApJ_ 707, 954–970.  10.1088/0004-637X/707/2/954
  * Pelupessy et al. (2006) Pelupessy, F. I., Papadopoulos, P. P., and van der Werf, P. (2006).  Incorporating the Molecular Gas Phase in Galaxy-sized Numerical Simulations: First Applications in Dwarf Galaxies.  _ApJ_ 645, 1024–1042.  10.1086/504366
  * Perets et al. (2007) Perets, H. B., Lederhendler, A., Biham, O., Vidali, G., Li, L., Swords, S., et al. (2007).  Molecular Hydrogen Formation on Amorphous Silicates under Interstellar Conditions.  _ApJL_ 661, L163–L166.  10.1086/518862
  * Peters et al. (2017) Peters, T., Naab, T., Walch, S., Glover, S. C. O., Girichidis, P., Pellegrini, E., et al. (2017).  The SILCC project - IV. Impact of dissociating and ionizing radiation on the interstellar medium and H$\alpha$ emission as a tracer of the star formation rate.  _MNRAS_ 466, 3293–3308.  10.1093/mnras/stw3216
  * Pettitt et al. (2015) Pettitt, A. R., Dobbs, C. L., Acreman, D. M., and Bate, M. R. (2015).  The morphology of the Milky Way - II. Reconstructing CO maps from disc galaxies with live stellar distributions.  _MNRAS_ 449, 3911–3926.  10.1093/mnras/stv600
  * Pettitt et al. (2014) Pettitt, A. R., Dobbs, C. L., Acreman, D. M., and Price, D. J. (2014).  The morphology of the Milky Way - I. Reconstructing CO maps from simulations in fixed potentials.  _MNRAS_ 444, 919–941.  10.1093/mnras/stu1075
  * Pettitt et al. (2020) Pettitt, A. R., Dobbs, C. L., Baba, J., Colombo, D., Duarte-Cabral, A., Egusa, F., et al. (2020).  How do different spiral arm models impact the ISM and GMC population?  _MNRAS_ 498, 1159–1174.  10.1093/mnras/staa2242
  * Pettitt et al. (2017) Pettitt, A. R., Tasker, E. J., Wadsley, J. W., Keller, B. W., and Benincasa, S. M. (2017).  Star formation and ISM morphology in tidally induced spiral structures.  _MNRAS_ 468, 4189–4204.  10.1093/mnras/stx736
  * Pfeffer et al. (2018) Pfeffer, J., Kruijssen, J. M. D., Crain, R. A., and Bastian, N. (2018).  The E-MOSAICS project: simulating the formation and co-evolution of galaxies and their star cluster populations.  _MNRAS_ 475, 4309–4346.  10.1093/mnras/stx3124
  * Phillips (1999) Phillips, J. P. (1999).  Rotation in molecular clouds.  _Astrophysics & Space Science_ 134, 241–254.  10.1051/aas:1999137
  * Pinto et al. (2012) Pinto, C., Verdini, A., Galli, D., and Velli, M. (2012).  Reflection and dissipation of Alfvén waves in interstellar clouds.  _A &A_ 544, A66.  10.1051/0004-6361/201219019
  * Piontek and Ostriker (2007) Piontek, R. A. and Ostriker, E. C. (2007).  Models of Vertically Stratified Two-Phase ISM Disks with MRI-Driven Turbulence.  _ApJ_ 663, 183–203.  10.1086/518103
  * Planck Collaboration et al. (2016) Planck Collaboration, Adam, R., Ade, P. A. R., Aghanim, N., Alves, M. I. R., Arnaud, M., et al. (2016).  Planck intermediate results. XXXII. The relative orientation between the magnetic field and structures traced by interstellar dust.  _A &A_ 586, A135.  10.1051/0004-6361/201425044
  * Rathjen et al. (2021) Rathjen, T.-E., Naab, T., Girichidis, P., Walch, S., Wünsch, R., Dinnbier, F., et al. (2021).  SILCC VI - Multiphase ISM structure, stellar clustering, and outflows with supernovae, stellar winds, ionizing radiation, and cosmic rays.  _MNRAS_ 504, 1039–1061.  10.1093/mnras/stab900
  * Reina-Campos et al. (2022) Reina-Campos, M., Keller, B. W., Kruijssen, J. M. D., Gensior, J., Trujillo-Gomez, S., Jeffreson, S. M. R., et al. (2022).  Introducing EMP-Pathfinder: modelling the simultaneous formation and evolution of stellar clusters in their host galaxies.  _MNRAS_ 517, 3144–3180.  10.1093/mnras/stac1934
  * Renaud et al. (2015) Renaud, F., Bournaud, F., Emsellem, E., Agertz, O., Athanassoula, E., Combes, F., et al. (2015).  Environmental regulation of cloud and star formation in galactic bars.  _MNRAS_ 454, 3299–3310.  10.1093/mnras/stv2223
  * Renaud et al. (2013) Renaud, F., Bournaud, F., Emsellem, E., Elmegreen, B., Teyssier, R., Alves, J., et al. (2013).  A sub-parsec resolution simulation of the Milky Way: global structure of the interstellar medium and properties of molecular clouds.  _MNRAS_ 436, 1836–1851.  10.1093/mnras/stt1698
  * Rice et al. (2016) Rice, T. S., Goodman, A. A., Bergin, E. A., Beaumont, C., and Dame, T. M. (2016).  A Uniform Catalog of Molecular Clouds in the Milky Way.  _ApJ_ 822, 52.  10.3847/0004-637X/822/1/52
  * Richings et al. (2014) Richings, A. J., Schaye, J., and Oppenheimer, B. D. (2014).  Non-equilibrium chemistry and cooling in the diffuse interstellar medium - II. Shielded gas.  _MNRAS_ 442, 2780–2796.  10.1093/mnras/stu1046
  * Rieder et al. (2022) Rieder, S., Dobbs, C., Bending, T., Liow, K. Y., and Wurster, J. (2022).  The formation and early evolution of embedded star clusters in spiral galaxies.  _MNRAS_ 509, 6155–6168.  10.1093/mnras/stab3425
  * Roberts and Hausman (1984) Roberts, J., W. W. and Hausman, M. A. (1984).  Spiral structure and star formation. I - Formation mechanisms and mean free paths.  _ApJ_ 277, 744–767.  10.1086/161746
  * Roberts and Stewart (1987) Roberts, J., William W. and Stewart, G. R. (1987).  The Role of Orbital Dynamics and Cloud-Cloud Collisions in the Formation of Giant Molecular Clouds in Global Spiral Structures.  _ApJ_ 314, 10.  10.1086/165035
  * Roberts (1969) Roberts, W. W. (1969).  Large-Scale Shock Formation in Spiral Galaxies and its Implications on Star Formation.  _ApJ_ 158, 123.  10.1086/150177
  * Robertson and Kravtsov (2008) Robertson, B. E. and Kravtsov, A. V. (2008).  Molecular Hydrogen and Global Star Formation Relations in Galaxies.  _ApJ_ 680, 1083–1111.  10.1086/587796
  * Rodriguez-Fernandez and Combes (2008) Rodriguez-Fernandez, N. J. and Combes, F. (2008).  Gas flow models in the Milky Way embedded bars.  _A &A_ 489, 115–133.  10.1051/0004-6361:200809644
  * Rogers and Pittard (2013) Rogers, H. and Pittard, J. M. (2013).  Feedback from winds and supernovae in massive stellar clusters - I. Hydrodynamics.  _MNRAS_ 431, 1337–1351.  10.1093/mnras/stt255
  * Roman-Duval et al. (2010) Roman-Duval, J., Jackson, J. M., Heyer, M., Rathborne, J., and Simon, R. (2010).  Physical Properties and Galactic Distribution of Molecular Clouds Identified in the Galactic Ring Survey.  _ApJ_ 723, 492–507.  10.1088/0004-637X/723/1/492
  * Rosdahl et al. (2015) Rosdahl, J., Schaye, J., Teyssier, R., and Agertz, O. (2015).  Galaxies that shine: radiation-hydrodynamical simulations of disc galaxies.  _MNRAS_ 451, 34–58.  10.1093/mnras/stv937
  * Rosen and Bregman (1995) Rosen, A. and Bregman, J. N. (1995).  Global Models of the Interstellar Medium in Disk Galaxies.  _ApJ_ 440, 634.  10.1086/175303
  * Rosolowsky (2007) Rosolowsky, E. (2007).  Giant Molecular Clouds in M31. I. Molecular Cloud Properties.  _ApJ_ 654, 240–251.  10.1086/509249
  * Rosolowsky et al. (2003) Rosolowsky, E., Engargiola, G., Plambeck, R., and Blitz, L. (2003).  Giant Molecular Clouds in M33. II. High-Resolution Observations.  _ApJ_ 599, 258–274.  10.1086/379166
  * Ruffle and Herbst (2000) Ruffle, D. P. and Herbst, E. (2000).  New models of interstellar gas-grain chemistry - I. Surface diffusion rates.  _MNRAS_ 319, 837–850.  10.1046/j.1365-8711.2000.03911.x
  * Salem and Bryan (2014) Salem, M. and Bryan, G. L. (2014).  Cosmic ray driven outflows in global galaxy disc models.  _MNRAS_ 437, 3312–3330.  10.1093/mnras/stt2121
  * Sales et al. (2014) Sales, L. V., Marinacci, F., Springel, V., and Petkova, M. (2014).  Stellar feedback by radiation pressure and photoionization.  _MNRAS_ 439, 2990–3006.  10.1093/mnras/stu155
  * Salo and Laurikainen (2000a) Salo, H. and Laurikainen, E. (2000a).  N-body model for M51 - I. Multiple encounter versus single passage?  _MNRAS_ 319, 377–392.  10.1046/j.1365-8711.2000.03650.x
  * Salo and Laurikainen (2000b) Salo, H. and Laurikainen, E. (2000b).  N-body model for M51 - II. Inner structure.  _MNRAS_ 319, 393–413.  10.1046/j.1365-8711.2000.03651.x
  * Sanders et al. (1984) Sanders, D. B., Solomon, P. M., and Scoville, N. Z. (1984).  Giant molecular clouds in the Galaxy. I - The axisymmetric distribution of H2.  _ApJ_ 276, 182–203.  10.1086/161602
  * Sano et al. (2018) Sano, H., Yamane, Y., Tokuda, K., Fujii, K., Tsuge, K., Nagaya, T., et al. (2018).  Molecular Clouds Associated with the Type Ia SNR N103B in the Large Magellanic Cloud.  _ApJ_ 867, 7.  10.3847/1538-4357/aae07c
  * Santillán et al. (2007) Santillán, A., Sánchez-Salcedo, F. J., and Franco, J. (2007).  Exploring Cloudy Gas Accretion as a Source of Interstellar Turbulence in the Outskirts of Disks.  _ApJL_ 662, L19–L22.  10.1086/519247
  * Sartorio et al. (2021) Sartorio, N. S., Vandenbroucke, B., Falceta-Goncalves, D., and Wood, K. (2021).  Photoionization feedback in turbulent molecular clouds.  _MNRAS_ 500, 1833–1843.  10.1093/mnras/staa3380
  * Schaye (2004) Schaye, J. (2004).  Star Formation Thresholds and Galaxy Edges: Why and Where.  _ApJ_ 609, 667–682.  10.1086/421232
  * Schmidt et al. (2013) Schmidt, W., Collins, D. C., and Kritsuk, A. G. (2013).  Local support against gravity in magnetoturbulent fluids.  _MNRAS_ 431, 3196–3215.  10.1093/mnras/stt399
  * Schneider et al. (2023) Schneider, N., Bonne, L., Bontemps, S., Kabanovic, S., Simon, R., Ossenkopf-Okada, V., et al. (2023).  Ionized carbon as a tracer of the assembly of interstellar clouds.  _Nature Astronomy_ 7, 546–556.  10.1038/s41550-023-01901-5
  * Seifried et al. (2020) Seifried, D., Haid, S., Walch, S., Borchert, E. M. A., and Bisbas, T. G. (2020).  SILCC-Zoom: H2 and CO-dark gas in molecular clouds - the impact of feedback and magnetic fields.  _MNRAS_ 492, 1465–1483.  10.1093/mnras/stz3563
  * Seifried et al. (2017) Seifried, D., Walch, S., Girichidis, P., Naab, T., Wünsch, R., Klessen, R. S., et al. (2017).  SILCC-Zoom: the dynamic and chemical evolution of molecular clouds.  _MNRAS_ 472, 4797–4818.  10.1093/mnras/stx2343
  * Seifried et al. (2018) Seifried, D., Walch, S., Haid, S., Girichidis, P., and Naab, T. (2018).  Is Molecular Cloud Turbulence Driven by External Supernova Explosions?  _ApJ_ 855, 81.  10.3847/1538-4357/aaacff
  * Sellwood (1980) Sellwood, J. A. (1980).  Galaxy models with live halos.  _A &A_ 89, 296–307
  * Sellwood and Carlberg (1984) Sellwood, J. A. and Carlberg, R. G. (1984).  Spiral instabilities provoked by accretion and star formation.  _ApJ_ 282, 61–74.  10.1086/162176
  * Semczuk et al. (2018) Semczuk, M., Łokas, E. L., Salomon, J.-B., Athanassoula, E., and D’Onghia, E. (2018).  Tidally Induced Morphology of M33 in Hydrodynamical Simulations of Its Recent Interaction with M31.  _ApJ_ 864, 34.  10.3847/1538-4357/aad4ae
  * Shetty et al. (2010) Shetty, R., Collins, D. C., Kauffmann, J., Goodman, A. A., Rosolowsky, E. W., and Norman, M. L. (2010).  The Effect of Projection on Derived Mass-Size and Linewidth-Size Relationships.  _ApJ_ 712, 1049–1056.  10.1088/0004-637X/712/2/1049
  * Shetty et al. (2011) Shetty, R., Glover, S. C., Dullemond, C. P., Ostriker, E. C., Harris, A. I., and Klessen, R. S. (2011).  Modelling CO emission - II. The physical characteristics that determine the X factor in Galactic molecular clouds.  _MNRAS_ 415, 3253–3274.  10.1111/j.1365-2966.2011.18937.x
  * Shetty et al. (2007) Shetty, R., Vogel, S. N., Ostriker, E. C., and Teuben, P. J. (2007).  Kinematics of Spiral-Arm Streaming in M51.  _ApJ_ 665, 1138–1158.  10.1086/520037
  * Shin et al. (2017) Shin, J., Kim, S. S., Baba, J., Saitoh, T. R., Hwang, J.-S., Chun, K., et al. (2017).  Hydrodynamic Simulations of the Central Molecular Zone with a Realistic Galactic Potential.  _ApJ_ 841, 74.  10.3847/1538-4357/aa7061
  * Shu et al. (1987) Shu, F. H., Adams, F. C., and Lizano, S. (1987).  Star formation in molecular clouds: observation and theory.  _ARA &A_ 25, 23–81.  10.1146/annurev.aa.25.090187.000323
  * Sillero et al. (2021) Sillero, E., Tissera, P. B., Lambas, D. G., Bovino, S., Schleicher, D. R., Grassi, T., et al. (2021).  Modelling H2 and its effects on star formation using a joint implementation of GADGET-3 and KROME.  _MNRAS_ 504, 2325–2345.  10.1093/mnras/stab1015
  * Sills et al. (2023) Sills, A., Rieder, S., Buckner, A. S. M., Hacar, A., Portegies Zwart, S., and Teixeira, P. S. (2023).  Using molecular gas observations to guide initial conditions for star cluster simulations.  _MNRAS_ 519, 4142–4151.  10.1093/mnras/stac3745
  * Skarbinski et al. (2023) Skarbinski, M., Jeffreson, S. M. R., and Goodman, A. A. (2023).  Building the molecular cloud population: the role of cloud mergers.  _MNRAS_ 519, 1887–1898.  10.1093/mnras/stac3627
  * Slyz et al. (2005) Slyz, A. D., Devriendt, J. E. G., Bryan, G., and Silk, J. (2005).  Towards simulating star formation in the interstellar medium.  _MNRAS_ 356, 737–752.  10.1111/j.1365-2966.2004.08494.x
  * Smilgys and Bonnell (2017) Smilgys, R. and Bonnell, I. A. (2017).  Formation of stellar clusters.  _MNRAS_ 472, 4982–4991.  10.1093/mnras/stx2396
  * Smith et al. (2022) Smith, J. D., Dale, J. E., Jaffa, S. E., and Krause, M. G. H. (2022).  Star cluster formation in clouds with externally driven turbulence.  _MNRAS_ 516, 4212–4219.  10.1093/mnras/stac2295
  * Smith et al. (2014) Smith, R. J., Glover, S. C. O., Clark, P. C., Klessen, R. S., and Springel, V. (2014).  CO-dark gas and molecular filaments in Milky Way-type galaxies.  _MNRAS_ 441, 1628–1645.  10.1093/mnras/stu616
  * Smith et al. (2023) Smith, R. J., Tress, R., Soler, J. D., Klessen, R. S., Glover, S. C. O., Hennebelle, P., et al. (2023).  On the distribution of the CNM in spiral galaxies.  _MNRAS_ 10.1093/mnras/stad1537
  * Smith et al. (2020) Smith, R. J., Treß, R. G., Sormani, M. C., Glover, S. C. O., Klessen, R. S., Clark, P. C., et al. (2020).  The Cloud Factory I: Generating resolved filamentary molecular clouds from galactic-scale forces.  _MNRAS_ 492, 1594–1613.  10.1093/mnras/stz3328
  * Smoluchowski (1981) Smoluchowski, R. (1981).  Rate of H2 formation on amorphous grains.  _Ap &SS_ 75, 353–363.  10.1007/BF00648648
  * Snow (1983) Snow, T. P. (1983).  On the possibly low H2 formation rate in dense clouds.  _ApJL_ 269, L57–L59.  10.1086/184055
  * Soler et al. (2013) Soler, J. D., Hennebelle, P., Martin, P. G., Miville-Deschênes, M. A., Netterfield, C. B., and Fissel, L. M. (2013).  An Imprint of Molecular Cloud Magnetization in the Morphology of the Dust Polarized Emission.  _ApJ_ 774, 128.  10.1088/0004-637X/774/2/128
  * Sormani et al. (2018) Sormani, M. C., Sobacchi, E., Fragkoudi, F., Ridley, M., Treß, R. G., Glover, S. C. O., et al. (2018).  A dynamical mechanism for the origin of nuclear rings.  _MNRAS_ 481, 2–19.  10.1093/mnras/sty2246
  * Sormani et al. (2019) Sormani, M. C., Treß, R. G., Glover, S. C. O., Klessen, R. S., Barnes, A. T., Battersby, C. D., et al. (2019).  The geometry of the gas surrounding the Central Molecular Zone: on the origin of localized molecular clouds with extreme velocity dispersions.  _MNRAS_ 488, 4663–4673.  10.1093/mnras/stz2054
  * Spilker et al. (2022) Spilker, A., Kainulainen, J., and Orkisz, J. (2022).  Bird’s eye view of molecular clouds in the Milky Way. II. Cloud kinematics from subparsec to kiloparsec scales.  _A &A_ 667, A110.  10.1051/0004-6361/202244392
  * Stinson et al. (2013) Stinson, G. S., Brook, C., Macciò, A. V., Wadsley, J., Quinn, T. R., and Couchman, H. M. P. (2013).  Making Galaxies In a Cosmological Context: the need for early stellar feedback.  _MNRAS_ 428, 129–140.  10.1093/mnras/sts028
  * Sun et al. (2020) Sun, J., Leroy, A. K., Schinnerer, E., Hughes, A., Rosolowsky, E., Querejeta, M., et al. (2020).  Molecular Gas Properties on Cloud Scales across the Local Star-forming Galaxy Population.  _ApJL_ 901, L8.  10.3847/2041-8213/abb3be
  * Sun et al. (2018) Sun, J., Leroy, A. K., Schruba, A., Rosolowsky, E., Hughes, A., Kruijssen, J. M. D., et al. (2018).  Cloud-scale Molecular Gas Properties in 15 Nearby Galaxies.  _ApJ_ 860, 172.  10.3847/1538-4357/aac326
  * Tanaka et al. (2007) Tanaka, K., Kamegai, K., Nagai, M., and Oka, T. (2007).  High-Resolution Mappings of the l$=$1 degree 3 Complex in Molecular Lines: Discovery of a Proto-Superbubble.  _PASJ_ 59, 323–333.  10.1093/pasj/59.2.323
  * Tanvir and Dale (2020) Tanvir, T. S. and Dale, J. E. (2020).  Collision between molecular clouds - I. The effect of the cloud virial ratio in head-on collisions.  _MNRAS_ 494, 246–258.  10.1093/mnras/staa665
  * Tasker (2011) Tasker, E. J. (2011).  Star Formation in Disk Galaxies. II. The Effect Of Star Formation and Photoelectric Heating on the Formation and Evolution of Giant Molecular Clouds.  _ApJ_ 730, 11–+.  10.1088/0004-637X/730/1/11
  * Tasker et al. (2008) Tasker, E. J., Bryan, G. L., and Tan, J. C. (2008).  Simulating the ISM in Global Disk Galaxies.  _ArXiv:0802.3944_
  * Tasker and Tan (2009) Tasker, E. J. and Tan, J. C. (2009).  Star Formation in Disk Galaxies. I. Formation and Evolution of Giant Molecular Clouds via Gravitational Instability and Cloud Collisions.  _ApJ_ 700, 358–375.  10.1088/0004-637X/700/1/358
  * Tenorio-Tagle and Palous (1987) Tenorio-Tagle, G. and Palous, J. (1987).  Giant-scale supernova remnants - The role of differential galactic rotation and the formation of molecular clouds.  _A &A_ 186, 287–294
  * Tielens and Hagen (1982) Tielens, A. G. G. M. and Hagen, W. (1982).  Model calculations of the molecular composition of interstellar grain mantles.  _A &A_ 114, 245–260
  * Tomisaka (1984) Tomisaka, K. (1984).  Coagulation of interstellar clouds in spiral gravitational potential and formation of giant molecular clouds.  _PASJ_ 36, 457–475
  * Toomre and Toomre (1972) Toomre, A. and Toomre, J. (1972).  Galactic Bridges and Tails.  _ApJ_ 178, 623–666.  10.1086/151823
  * Tress et al. (2020a) Tress, R. G., Smith, R. J., Sormani, M. C., Glover, S. C. O., Klessen, R. S., Mac Low, M.-M., et al. (2020a).  Simulations of the star-forming molecular gas in an interacting M51-like galaxy.  _MNRAS_ 492, 2973–2995.  10.1093/mnras/stz3600
  * Tress et al. (2020b) Tress, R. G., Smith, R. J., Sormani, M. C., Glover, S. C. O., Klessen, R. S., Mac Low, M.-M., et al. (2020b).  Simulations of the star-forming molecular gas in an interacting M51-like galaxy.  _MNRAS_ 492, 2973–2995.  10.1093/mnras/stz3600
  * Valentini et al. (2023) Valentini, M., Dolag, K., Borgani, S., Murante, G., Maio, U., Tornatore, L., et al. (2023).  Impact of H2-driven star formation and stellar feedback from low-enrichment environments on the formation of spiral galaxies.  _MNRAS_ 518, 1128–1147.  10.1093/mnras/stac2110
  * Van Loo et al. (2013) Van Loo, S., Butler, M. J., and Tan, J. C. (2013).  Kiloparsec-scale Simulations of Star Formation in Disk Galaxies. I. The Unmagnetized and Zero-feedback Limit.  _ApJ_ 764, 36.  10.1088/0004-637X/764/1/36
  * Vandenbroucke and Wood (2019) Vandenbroucke, B. and Wood, K. (2019).  Radiation hydrodynamics simulations of the evolution of the diffuse ionized gas in disc galaxies.  _MNRAS_ 488, 1977–1986.  10.1093/mnras/stz1841
  * Vázquez-Semadeni et al. (2011) Vázquez-Semadeni, E., Banerjee, R., Gómez, G. C., Hennebelle, P., Duffin, D., and Klessen, R. S. (2011).  Molecular cloud evolution - IV. Magnetic fields, ambipolar diffusion and the star formation efficiency.  _MNRAS_ 414, 2511–2527.  10.1111/j.1365-2966.2011.18569.x
  * Vázquez-Semadeni et al. (2000) Vázquez-Semadeni, E., Gazol, A., and Scalo, J. (2000).  Is Thermal Instability Significant in Turbulent Galactic Gas?  _ApJ_ 540, 271–285.  10.1086/309318
  * Wada and Koda (2004) Wada, K. and Koda, J. (2004).  Instabilities of spiral shocks - I. Onset of wiggle instability and its mechanism.  _MNRAS_ 349, 270–280
  * Wada et al. (2002) Wada, K., Meurer, G., and Norman, C. A. (2002).  Gravity-driven Turbulence in Galactic Disks.  _ApJ_ 577, 197–205.  10.1086/342151
  * Wada and Norman (1999) Wada, K. and Norman, C. A. (1999).  The Global Structure and Evolution of a Self-Gravitating Multiphase Interstellar Medium in a Galactic Disk.  _ApJL_ 516, L13–L16.  10.1086/311987
  * Wada and Norman (2001) Wada, K. and Norman, C. A. (2001).  Numerical Models of the Multiphase Interstellar Matter with Stellar Energy Feedback on a Galactic Scale.  _ApJ_ 547, 172–186.  10.1086/318344
  * Wada et al. (2000) Wada, K., Spaans, M., and Kim, S. (2000).  Formation of Cavities, Filaments, and Clumps by the Nonlinear Development of Thermal and Gravitational Instabilities in the Interstellar Medium under Stellar Feedback.  _ApJ_ 540, 797–807.  10.1086/309347
  * Wada et al. (1994) Wada, K., Taniguchi, Y., Habe, A., and Hasegawa, T. (1994).  Gasdynamical Approach to a Face-on View of the Milky Way.  _ApJL_ 437, L123.  10.1086/187698
  * Walch et al. (2015) Walch, S., Girichidis, P., Naab, T., Gatto, A., Glover, S. C. O., Wünsch, R., et al. (2015).  The SILCC (SImulating the LifeCycle of molecular Clouds) project - I. Chemical evolution of the supernova-driven ISM.  _MNRAS_ 454, 238–268.  10.1093/mnras/stv1975
  * Wareing et al. (2016) Wareing, C. J., Pittard, J. M., Falle, S. A. E. G., and Van Loo, S. (2016).  Magnetohydrodynamical simulation of the formation of clumps and filaments in quiescent diffuse medium by thermal instability.  _MNRAS_ 459, 1803–1818.  10.1093/mnras/stw581
  * Watkins et al. (2023a) Watkins, E. J., Barnes, A. T., Henny, K., Kim, H., Kreckel, K., Meidt, S. E., et al. (2023a).  PHANGS-JWST First Results: A Statistical View on Bubble Evolution in NGC 628.  _ApJL_ 944, L24.  10.3847/2041-8213/aca6e4
  * Watkins et al. (2023b) Watkins, E. J., Kreckel, K., Groves, B., Glover, S. C. O., Whitmore, B. C., Leroy, A. K., et al. (2023b).  Quantifying the energetics of molecular superbubbles in PHANGS galaxies.  _A &A_ 676, A67.  10.1051/0004-6361/202346075
  * Whitworth et al. (2023) Whitworth, D. J., Smith, R. J., Klessen, R. S., Mac Low, M.-M., Glover, S. C. O., Tress, R., et al. (2023).  Magnetic fields do not suppress global star formation in low metallicity dwarf galaxies.  _MNRAS_ 520, 89–106.  10.1093/mnras/stad105
  * Wibking and Krumholz (2023) Wibking, B. D. and Krumholz, M. R. (2023).  The global structure of magnetic fields and gas in simulated Milky Way-analogue galaxies.  _MNRAS_ 521, 5972–5990.  10.1093/mnras/stac2648
  * Wiener et al. (2013) Wiener, J., Zweibel, E. G., and Oh, S. P. (2013).  Cosmic Ray Heating of the Warm Ionized Medium.  _ApJ_ 767, 87.  10.1088/0004-637X/767/1/87
  * Williamson et al. (2014) Williamson, D. J., Thacker, R. J., Wurster, J., and Gibson, B. K. (2014).  Cloud angular momentum and effective viscosity in global SPH simulations with feedback.  _MNRAS_ 442, 3674–3685.  10.1093/mnras/stu1121
  * Wu et al. (2017) Wu, B., Tan, J. C., Christie, D., Nakamura, F., Van Loo, S., and Collins, D. (2017).  GMC Collisions as Triggers of Star Formation. III. Density and Magnetically Regulated Star Formation.  _ApJ_ 841, 88.  10.3847/1538-4357/aa6ffa
  * Zucker et al. (2022) Zucker, C., Goodman, A. A., Alves, J., Bialy, S., Foley, M., Speagle, J. S., et al. (2022).  Star formation near the Sun is driven by expansion of the Local Bubble.  _Nature_ 601, 334–337.  10.1038/s41586-021-04286-5


# Novel Local Characteristic Decomposition Based Path-Conservative Central-
Upwind Schemes

Shaoshuai Chu, Michael Herty, and Alexander Kurganov Department of Mathematics
and Shenzhen International Center for Mathematics, Southern University of
Science and Technology, Shenzhen, 518055, China;
<EMAIL_ADDRESS>of Mathematics, RWTH Aachen University,
52056 Aachen, Germany<EMAIL_ADDRESS>of Mathematics,
Shenzhen International Center for Mathematics, and Guangdong Provincial Key
Laboratory of Computational Science and Material Design, Southern University
of Science and Technology, Shenzhen, 518055, China<EMAIL_ADDRESS>
###### Abstract

We introduce local characteristic decomposition based path-conservative
central-upwind schemes for (nonconservative) hyperbolic systems of balance
laws. The proposed schemes are made to be well-balanced via a flux
globalization approach, in which source terms are incorporated into the
fluxes: This helps to enforce the well-balanced property when the resulting
quasi-conservative system is solved using the local characteristic
decomposition based central-upwind scheme recently introduced in [A. Chertock,
S. Chu, M. Herty, A. Kurganov, and M. Lukáčová-Medvid’ová, J. Comput. Phys.,
473 (2023), Paper No. 111718]. Nonconservative product terms are also
incorporated into the global fluxes using a path-conservative technique. We
illustrate the performance of the developed schemes by applying them to one-
and two-dimensional compressible multifluid systems and thermal rotating
shallow water equations.

Key words: Local characteristic decomposition, path-conservative central-
upwind schemes, flux globalization, compressible multifluids, thermal rotating
shallow water equations.

AMS subject classification: 76M12, 65M08, 35L65, 86-08, 76T99.

##  1 Introduction

This paper is focused on the development of highly accurate and robust
numerical methods for the nonconservative hyperbolic systems of balance laws,
which in the one-dimensional (1-D) case read as

$\bm{U}_{t}+\bm{F}(\bm{U})_{x}=B(\bm{U})\bm{U}_{x}+\bm{S}(\bm{U}),$ (1.1)

where $x$ is a spatial variable, $t$ is the time, $\bm{U}\in\mathbb{R}^{d}$ is
a vector of unknown functions, $\bm{F}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ are
nonlinear fluxes, $B(\bm{U})\in\mathbb{R}^{d\times d}$, and
$\bm{S}(\bm{U}):\mathbb{R}^{d}\to\mathbb{R}^{d}$ are source terms.

Designing a numerical scheme for (1.1) is a challenging task for a number of
reasons. First, it is well-known that solutions of (1.1) may develop
complicated nonsmooth structures even for smooth initial data. Second, the
presence of nonconservative products (when $B\neq 0$) brings many challenges
for studying these systems both theoretically and numerically. In general,
when $\bm{U}$ is discontinuous, weak solutions of (1.1) cannot be understood
in the sense of distributions. Instead, one can introduce Borel measure
solutions [19, 30, 31]. This concept of weak solutions was utilized to develop
path-conservative finite-volume (FV) methods; see, e.g., [35, 6, 7, 8, 36, 38,
1, 9] and reference therein. Third, a novel scheme should be well-balanced
(WB) in the sense that it should be capable of exactly preserving some of the
physically relevant steady-state solutions satisfying

$\bm{F}(\bm{U})_{x}-B(\bm{U})\bm{U}_{x}-\bm{S}(\bm{U})=M(\bm{U})\bm{E}(\bm{U})_{x}=\bm{0}.$

Here, $M\in\mathbb{R}^{d\times d}$ and $\bm{E}$ is the vector of equilibrium
variables, which are constant at the steady states provided $M$ is invertible;
see [26].

In this paper, we focus on one of the popular classes of FV methods, the
Riemann-problem-solver-free central-upwind (CU) schemes, which were originally
introduced in [28, 29, 25] as a “black-box” solver for general
multidimensional hyperbolic systems of conservation laws and then successfully
applied to a variety of hyperbolic systems. In [8], the CU schemes were first
reformulated in a path-conservative framework and then extended to the 1-D
nonconservative systems (1.1), resulting in the second-order path-conservative
CU (PCCU) schemes. These schemes are WB in the sense that they are capable of
preserving simple steady-state solutions like “lake-at-rest” of several
shallow water models.

In order to design PCCU schemes capable of preserving more complicated steady
states, path-conservative techniques were incorporated into the flux
globalization framework: This results in a new class of flux globalization
based WB PCCU schemes introduced in [26]. A flux globalization approach, which
was introduced in [14, 5, 20, 21, 33], relies on the following quasi-
conservative form of (1.1):

$\bm{U}_{t}+\bm{K}(\bm{U})_{x}=\bm{0},$ (1.2)

where $\bm{K}(\bm{U})$ is a global flux

$\bm{K}(\bm{U})=\bm{F}(\bm{U})-\bm{W}(\bm{U}),\quad\bm{W}(\bm{U})=\int\limits^{x}_{\widehat{x}}B(\bm{U}){\bm{U}}_{\xi}(\xi,t)+S(\bm{U}(\xi,t))\,{\rm
d}\xi,$ (1.3)

and ${\widehat{x}}$ is an arbitrary number. The flux globalization based WB
PCCU schemes are obtained by applying the CU numerical flux to the semi-
discretization of (1.2) and by using the path-conservative technique for the
evaluation of the integrals in (1.3). These schemes have been applied to a
variety of hyperbolic systems in [2, 3, 4, 10, 26, 18].

Even though the flux globalization based WB PCCU schemes from [26] are
accurate, efficient, and robust, higher resolution of the numerical solutions
can be achieved by reducing numerical dissipation present in those schemes.
This can be done, for example, by using the local characteristic decomposition
(LCD) based CU (LCD-CU) numerical fluxes instead of the original CU ones
employed in [2, 3, 4, 10, 26]. The LCD-CU schemes have been recently
introduced in [12] for hyperbolic systems of conservation laws and we use the
same LCD based framework to enhance the performance of the flux globalization
based WB PCCU schemes. The proposed scheme combines for the first time the
modifications enhancing the performance of the CU schemes for both
conservative and nonconservative hyperbolic systems in multiple space
dimensions. We derive the new LCD-PCCU schemes and apply them to two
particular systems. First, we consider the $\gamma$-based compressible
multifluid systems, which can be, in principle, written in the conservative
form, but we use its nonconservative formulation to achieve high resolution of
material interfaces. Second, we consider a challenging thermal rotating
shallow water (TRSW) model, which contains geometric source term due to a
nonflat bottom topography, source terms due to Coriolis forces, and
generically nonconservative momentum exchange terms.

The rest of the paper is organized as follows. In §2, we briefly describe the
1-D second-order LCD-CU scheme. In §3, we propose the novel scheme by
considering the nonconservative system (1.1) and design a 1-D flux
globalization based WB LCD-PCCU scheme. In §4, we extend the 1-D LCD-PCCU
scheme to the two-dimensional (2-D) case. Finally, in §5, we apply the
developed LCD-PCCU schemes to the $\gamma$-based compressible multifluid
systems and TRSW equations. We illustrate the performance of the scheme on a
number of numerical examples.

##  2 1-D LCD-CU Scheme

In this section, we briefly describe the second-order semi-discrete LCD-CU
scheme, which is a slightly modified version of the scheme from [12]. For the
difference between the schemes, we refer the reader to Remarks 2.1–2.3 below.

We consider a 1-D hyperbolic system of conservative laws

$\bm{U}_{t}+\bm{F}(\bm{U})_{x}=\bm{0},$ (2.1)

cover the computational domain with the uniform cells
$C_{j}:=[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]$ of size $\Delta x$ centered at
$x_{j}=(x_{j-\frac{1}{2}}+x_{j+\frac{1}{2}})/2$, and denote by
$\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j}(t)$ the computed cell
averages of $\bm{U}(\cdot,t)$,

$\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j}(t):\approx\frac{1}{\Delta x}\int\limits_{C_{j}}\bm{U}(x,t)\,{\rm
d}x,\qquad j=1,\ldots,N,$

which are assumed to be available at a certain time $t\geq 0$. Note that both
$\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j}(t)$ and many other
indexed quantities below are time-dependent, but from here on, we suppress
their time-dependence for the sake of brevity.

The computed cell averages (2.1) are evolved in time by numerically solving
the following system of ODEs:

$\frac{{\rm d}\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j}}{{\rm
d}t}=-\frac{\bm{\mathcal{F}}^{\rm LCD}_{j+\frac{1}{2}}-\bm{\mathcal{F}}^{\rm
LCD}_{j-\frac{1}{2}}}{\Delta x},$

where the LCD-CU numerical fluxes are given by

$\bm{\mathcal{F}}^{\rm LCD}_{j+\frac{1}{2}}=R_{j+\frac{1}{2}}P^{\rm
LCD}_{j+\frac{1}{2}}R^{-1}_{j+\frac{1}{2}}\bm{F}^{-}_{j+\frac{1}{2}}+R_{j+\frac{1}{2}}M^{\rm
LCD}_{j+\frac{1}{2}}R^{-1}_{j+\frac{1}{2}}\bm{F}^{+}_{j+\frac{1}{2}}+R_{j+\frac{1}{2}}Q^{\rm
LCD}_{j+\frac{1}{2}}R^{-1}_{j+\frac{1}{2}}\big{(}\bm{U}^{+}_{j+\frac{1}{2}}-\bm{U}^{-}_{j+\frac{1}{2}}\big{)}.$
(2.2)

Here, $\bm{F}^{\pm}_{j}=\bm{F}(\bm{U}^{\pm}_{j})$ and
$\bm{U}^{\pm}_{j+\frac{1}{2}}$ are the one-sided values of a piecewise linear
interpolant

$\widetilde{\bm{U}}(x)=\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j}+(\bm{U}_{x})_{j}(x-x_{j}),\quad x\in C_{j},$ (2.3)

at the cell interface $x=x_{j+\frac{1}{2}}$, namely

$\bm{U}^{-}_{j+\frac{1}{2}}=\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j}+\frac{\Delta
x}{2}(\bm{U}_{x})_{j},\quad\bm{U}^{+}_{j+\frac{1}{2}}=\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j+1}-\frac{\Delta x}{2}(\bm{U}_{x})_{j+1}.$ (2.4)

In order to ensure the non-oscillatory nature of this reconstruction, we
compute the slopes $(\bm{U}_{x})_{j}$ in (2.4) using a generalized minmod
limiter [32, 34, 40]:

$(\bm{U}_{x})_{j}={\rm minmod}\left(\theta\frac{\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j}-\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j-1}}{\Delta
x},\,\frac{\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j+1}-\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j-1}}{2\Delta
x},\,\theta\frac{\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j+1}-\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j}}{\Delta
x}\right),\quad\theta\in[1,2],$ (2.5)

applied in the component-wise manner. Here, the minmod function is defined as

${\rm
minmod}(z_{1},z_{2},\ldots):=\begin{cases}\min_{j}\\{z_{j}\\}&\mbox{if}~{}z_{j}>0\quad\forall\,j,\\\
\max_{j}\\{z_{j}\\}&\mbox{if}~{}z_{j}<0\quad\forall\,j,\\\
0&\mbox{otherwise},\end{cases}$ (2.6)

and the parameter $\theta$ can be used to control the amount of numerical
viscosity present in the resulting scheme. We have taken $\theta=1.3$ in all
of the numerical experiments reported in §5.

The diagonal matrices $P^{\rm LCD}_{j+\frac{1}{2}}$, $M^{\rm
LCD}_{j+\frac{1}{2}}$, and $Q^{\rm LCD}_{j+\frac{1}{2}}$ are given by

$\displaystyle P^{\rm LCD}_{j+\frac{1}{2}}={\rm diag}\big{(}(P^{\rm
LCD}_{1})_{j+\frac{1}{2}},\ldots,(P^{\rm
LCD}_{d})_{j+\frac{1}{2}}\big{)},\quad M^{\rm LCD}_{j+\frac{1}{2}}={\rm
diag}\big{(}(M^{\rm LCD}_{1})_{j+\frac{1}{2}},\ldots,(M^{\rm
LCD}_{d})_{j+\frac{1}{2}}\big{)},$ (2.7) $\displaystyle Q^{\rm
LCD}_{j+\frac{1}{2}}={\rm diag}\big{(}(Q^{\rm
LCD}_{1})_{j+\frac{1}{2}},\ldots,(Q^{\rm LCD}_{d})_{j+\frac{1}{2}}\big{)}$

with

$\displaystyle\big{(}(P^{\rm LCD}_{i})_{j+\frac{1}{2}},(M^{\rm
LCD}_{i})_{j+\frac{1}{2}},(Q^{\rm LCD}_{i})_{j+\frac{1}{2}}\big{)}$ (2.8)
$\displaystyle=\left\\{\begin{aligned}
&\frac{1}{(\lambda^{+}_{i})_{j+\frac{1}{2}}-(\lambda^{-}_{i})_{j+\frac{1}{2}}}\Big{(}(\lambda^{+}_{i})_{j+\frac{1}{2}},-(\lambda^{-}_{i})_{j+\frac{1}{2}},(\lambda^{+}_{i})_{j+\frac{1}{2}}(\lambda^{-}_{i})_{j+\frac{1}{2}}\Big{)}&&\mbox{if}~{}(\lambda^{+}_{i})_{j+\frac{1}{2}}-(\lambda^{-}_{i})_{j+\frac{1}{2}}>\varepsilon_{0},\\\\[1.29167pt]
&\frac{1}{a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}}\Big{(}a^{+}_{j+\frac{1}{2}},-a^{-}_{j+\frac{1}{2}},a^{+}_{j+\frac{1}{2}}a^{-}_{j+\frac{1}{2}}\Big{)}&&\mbox{if}~{}(\lambda^{+}_{i})_{j+\frac{1}{2}}-(\lambda^{-}_{i})_{j+\frac{1}{2}}\leq\varepsilon_{0}~{}~{}\mbox{and}~{}~{}a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}>\varepsilon_{0},\\\\[1.29167pt]
&\Big{(}{\frac{1}{2}},{\frac{1}{2}},0\Big{)}&&\mbox{if}~{}a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}\leq\varepsilon_{0},\end{aligned}\right.$

where

$\begin{aligned}
&(\lambda^{+}_{i})_{j+\frac{1}{2}}=\max\left\\{\lambda_{i}\big{(}A(\bm{U}^{-}_{j+\frac{1}{2}})\big{)},\,\lambda_{i}\big{(}A(\bm{U}^{+}_{j+\frac{1}{2}})\big{)},\,0\right\\},\\\
&(\lambda^{-}_{i})_{j+\frac{1}{2}}=\min\left\\{\lambda_{i}\big{(}A(\bm{U}^{-}_{j+\frac{1}{2}})\big{)},\,\lambda_{i}\big{(}A(\bm{U}^{+}_{j+\frac{1}{2}})\big{)},\,0\right\\},\end{aligned}\quad
i=1,\ldots,d.$

Here, $A(\bm{U})=\frac{\partial\bm{F}}{\partial\bm{U}}(\bm{U})$ is the
Jacobian and the matrices $R_{j+\frac{1}{2}}$ and $R^{-1}_{j+\frac{1}{2}}$ are
such that $R^{-1}_{j+\frac{1}{2}}\widehat{A}_{j+\frac{1}{2}}R_{j+\frac{1}{2}}$
is diagonal, where
$\widehat{A}_{j+\frac{1}{2}}=A(\widehat{\bm{U}}_{j+\frac{1}{2}})$ and
$\widehat{\bm{U}}_{j+\frac{1}{2}}$ is either a simple average
$(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j}+\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j+1})/2$ or another type of average of the $\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j}$
and $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j+1}$ states. The one-
sided local speeds of propagation $a^{\pm}_{j+\frac{1}{2}}$ can be estimated
using the largest and the smallest eigenvalues of $A$, for example, by taking

$\displaystyle
a^{+}_{j+\frac{1}{2}}=\max\left\\{\lambda_{d}\big{(}A(\bm{U}^{+}_{j+\frac{1}{2}})\big{)},\lambda_{d}\big{(}A(\bm{U}^{-}_{j+\frac{1}{2}})\big{)},0\right\\},\quad
a^{-}_{j+\frac{1}{2}}=\min\left\\{\lambda_{1}\big{(}A(\bm{U}^{+}_{j+\frac{1}{2}})\big{)},\lambda_{1}\big{(}A(\bm{U}^{-}_{j+\frac{1}{2}})\big{)},0\right\\}.$

Finally, $\varepsilon_{0}$ in (2.8) is a very small desingularization
constant, which is taken $\varepsilon_{0}=10^{-18}$ in all of the numerical
examples reported in §5.

###### Remark 2.1

The numerical flux (2.2) is identical to the LCD-CU numerical flux in [12],
which we have simplified using the fact that

$P^{\rm LCD}_{j+\frac{1}{2}}+M^{\rm LCD}_{j+\frac{1}{2}}=I.$ (2.9)

###### Remark 2.2

We perform a piecewise linear reconstruction for the conservative variables
$\bm{U}$ rather than the local characteristic variables as it was done in
[12]. As no substantial spurious oscillations have been observed in any of the
conducted numerical examples, there is no need to switch to the local
characteristic variables in this case. This allows the developed scheme to be
more efficient than the original LCD-CU scheme from [12].

###### Remark 2.3

The desingularization in (2.8) is performed differently than in [12]. Our
numerical experiments indicate that the current desingularization is more
robust compared with the previous desingularization strategy.

##  3 1-D Flux Globalization Based LCD-PCCU Scheme

In this section, we propose an extension of the LCD-CU scheme from §2 to the
nonconservative system (1.1) and design a 1-D flux globalization based LCD-
PCCU scheme.

We apply the LCD-CU scheme to the quasi-conservative system (1.2)–(1.3) to
obtain

$\frac{{\rm d}\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j}}{{\rm
d}t}=-\frac{\bm{\mathcal{K}}^{\rm LCD}_{j+\frac{1}{2}}-\bm{\mathcal{K}}^{\rm
LCD}_{j-\frac{1}{2}}}{\Delta x},$ (3.1)

where the numerical fluxes $\bm{\mathcal{K}}^{\rm LCD}_{j+\frac{1}{2}}$ are
given by (compare with (2.2))

$\bm{\mathcal{K}}^{\rm LCD}_{j+\frac{1}{2}}=R_{j+\frac{1}{2}}P^{\rm
LCD}_{j+\frac{1}{2}}R^{-1}_{j+\frac{1}{2}}\bm{K}^{-}_{j+\frac{1}{2}}+R_{j+\frac{1}{2}}M^{\rm
LCD}_{j+\frac{1}{2}}R^{-1}_{j+\frac{1}{2}}\bm{K}^{+}_{j+\frac{1}{2}}+R_{j+\frac{1}{2}}Q^{\rm
LCD}_{j+\frac{1}{2}}R^{-}_{j+\frac{1}{2}}\big{(}\breve{\bm{U}}^{+}_{j+\frac{1}{2}}-\breve{\bm{U}}^{-}_{j+\frac{1}{2}}\big{)}.$
(3.2)

Here, $\bm{K}^{\pm}_{j+\frac{1}{2}}$ are the global fluxes, which are computed
as the same as in [3, §2].

The global fluxes $\bm{K}^{\pm}_{j+\frac{1}{2}}$ are computed using the point
values $\bm{U}^{\pm}_{j+\frac{1}{2}}$, which, in order to ensure the WB
property of the resulting scheme, are to be obtained with the help of a
piecewise linear reconstruction (2.3)–(2.6) applied to the equilibrium
variables $\bm{E}$ (instead of $\bm{U}$):

$\bm{E}^{-}_{j+\frac{1}{2}}=\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{E}$\kern 0.0pt}}}_{j}+\frac{\Delta
x}{2}(\bm{E}_{x})_{j},\quad\bm{E}^{+}_{j+\frac{1}{2}}=\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{E}$\kern
0.0pt}}}_{j+1}-\frac{\Delta x}{2}(\bm{E}_{x})_{j+1},$

where $\bm{E}_{j}:=\bm{E}(\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j})$. We then obtain the
corresponding values $\bm{U}^{\pm}_{j+\frac{1}{2}}$ by solving the (nonlinear
systems of) equations

$\bm{E}(\bm{U}^{+}_{j+\frac{1}{2}})=\bm{E}^{+}_{j+\frac{1}{2}}\quad\mbox{and}\quad\bm{E}(\bm{U}^{-}_{j+\frac{1}{2}})=\bm{E}^{-}_{j+\frac{1}{2}}$
(3.3)

for $\bm{U}^{+}_{j+\frac{1}{2}}$ and $\bm{U}^{-}_{j+\frac{1}{2}}$,
respectively. We emphasize that the values $\bm{U}^{\pm}_{j+\frac{1}{2}}$
cannot be used in (3.2) instead of $\breve{\bm{U}}^{\pm}_{j+\frac{1}{2}}$ as
$\bm{U}^{+}_{j+\frac{1}{2}}$ may not be equal to $\bm{U}^{-}_{j+\frac{1}{2}}$
at discrete steady states. We therefore use another set of point values
$\breve{\bm{U}}^{\pm}_{j+\frac{1}{2}}$, which are obtained by solving modified
versions of (3.3):

$\breve{\bm{E}}(\breve{\bm{U}}^{+}_{j+\frac{1}{2}})=\bm{E}^{+}_{j+\frac{1}{2}}\quad\mbox{and}\quad\breve{\bm{E}}(\breve{\bm{U}}^{-}_{j+\frac{1}{2}})=\bm{E}^{-}_{j+\frac{1}{2}}.$

Here, the functions $\breve{\bm{E}}$ are the modifications of the functions
$\bm{E}$, which are made in such a way that
$\breve{\bm{U}}^{+}_{j+\frac{1}{2}}=\breve{\bm{U}}^{-}_{j+\frac{1}{2}}$ as
long as $\bm{E}^{+}_{j+\frac{1}{2}}=\bm{E}^{-}_{j+\frac{1}{2}}$ and, at the
same time,
$\breve{\bm{U}}^{+}_{j+\frac{1}{2}}=\bm{U}^{+}_{j+\frac{1}{2}}+{\cal
O}((\Delta x)^{2})$ and
$\breve{\bm{U}}^{-}_{j+\frac{1}{2}}=\bm{U}^{-}_{j+\frac{1}{2}}+{\cal
O}((\Delta x)^{2})$ for smooth solutions; see [10, 3, 2, 4, 26] for examples
of such modifications.

The diagonal matrices $P^{\rm LCD}_{j+\frac{1}{2}}$, $M^{\rm
LCD}_{j+\frac{1}{2}}$, and $Q^{\rm LCD}_{j+\frac{1}{2}}$ in (3.2) are defined
by (2.7)–(2.8), but with

$\begin{aligned}
(\lambda^{+}_{i})_{j+\frac{1}{2}}&=\max\left\\{\lambda_{i}\big{(}{\cal
A}(\bm{U}^{-}_{j+\frac{1}{2}})\big{)},\,\lambda_{i}\big{(}{\cal
A}(\bm{U}^{+}_{j+\frac{1}{2}})\big{)},\,0\right\\},\\\
(\lambda^{-}_{i})_{j+\frac{1}{2}}&=\min\left\\{\lambda_{i}\big{(}{\cal
A}(\bm{U}^{-}_{j+\frac{1}{2}})\big{)},\,\lambda_{i}\big{(}{\cal
A}(\bm{U}^{+}_{j+\frac{1}{2}})\big{)},\,0\right\\},\end{aligned}\quad
i=1,\ldots,d,$

where ${\cal
A}(\bm{U})=\frac{\partial\bm{F}}{\partial\bm{U}}(\bm{U})-B(\bm{U})$ and the
matrices $R_{j+\frac{1}{2}}$ and $R^{-1}_{j+\frac{1}{2}}$ are such that
$R^{-1}_{j+\frac{1}{2}}\widehat{\cal A}_{j+\frac{1}{2}}R_{j+\frac{1}{2}}$ is
diagonal for $\widehat{\cal A}_{j+\frac{1}{2}}={\cal
A}(\widehat{\bm{U}}_{j+\frac{1}{2}})$. Finally, the one-sided local speeds of
propagation $a^{\pm}_{j+\frac{1}{2}}$ are estimated by

$\displaystyle a^{+}_{j+\frac{1}{2}}=\max\left\\{\lambda_{d}\big{(}{\cal
A}(\bm{U}^{+}_{j+\frac{1}{2}})\big{)},\lambda_{d}\big{(}{\cal
A}(\bm{U}^{-}_{j+\frac{1}{2}})\big{)},\,0\right\\},\quad
a^{-}_{j+\frac{1}{2}}=\min\left\\{\lambda_{1}\big{(}{\cal
A}(\bm{U}^{+}_{j+\frac{1}{2}})\big{)},\lambda_{1}\big{(}{\cal
A}(\bm{U}^{-}_{j+\frac{1}{2}})\big{)},\,0\right\\}.$

###### Remark 3.1

The introduced 1-D semi-discrete flux globalization based LCD-PCCU scheme is
WB. In order to prove this, we only need to show that if the discrete data are
at a steady state, namely, if $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern
0.0pt}}}_{j}\equiv\widehat{E}=\texttt{Const}$, then the RHS of (3.1) vanishes.
Following the proof of [3, Theorem 2.1], we obtain

$\bm{K}^{+}_{j+\frac{1}{2}}=\bm{K}^{-}_{j+\frac{1}{2}}=\bm{K}^{+}_{j-\frac{1}{2}}=\widehat{\bm{K}}=\texttt{Const},\quad\forall
j.$ (3.4)

Substituting (3.4) and (2.9) into (3.2) and taking into account the fact that
at the steady states
$\breve{\bm{U}}^{+}_{j+\frac{1}{2}}=\breve{\bm{U}}^{-}_{j+\frac{1}{2}},~{}\forall
j$, yields

$\bm{\mathcal{K}}^{\rm
LCD}_{j+\frac{1}{2}}=\widehat{\bm{K}}+R_{j+\frac{1}{2}}Q^{\rm
LCD}_{j+\frac{1}{2}}R^{-}_{j+\frac{1}{2}}\big{(}\breve{\bm{U}}^{+}_{j+\frac{1}{2}}-\breve{\bm{U}}^{-}_{j+\frac{1}{2}}\big{)}=\widehat{\bm{K}},\quad\forall
j,$

which ensures that the RHS of (3.1) vanishes.

##  4 2-D Flux Globalization Based LCD-PCCU Scheme

In this section, we propose a novel 2-D LCD-PCCU scheme for

$\bm{U}_{t}+\bm{F}(\bm{U})_{x}+\bm{G}(\bm{U})_{y}=B(\bm{U})\bm{U}_{x}+C(\bm{U})\bm{U}_{y}+\bm{S}^{x}(\bm{U})+\bm{S}^{y}(\bm{U}),$
(4.1)

where $\bm{G}:\mathbb{R}^{d}\to\mathbb{R}^{d}$,
$C(\bm{U})\in\mathbb{R}^{d\times d}$,
$\bm{S}^{x}(\bm{U}):\mathbb{R}^{d}\to\mathbb{R}^{d}$, and
$\bm{S}^{y}(\bm{U}):\mathbb{R}^{d}\to\mathbb{R}^{d}$. To this end, we design
the 2-D scheme in a “dimension-by-dimension” manner.

As in the 1-D case, we first rewrite the general 2-D nonconservative system
(4.1) in the following quasi-conservative form:

$\bm{U}_{t}+\bm{K}(\bm{U})_{x}+\bm{L}(\bm{U})_{y}=\bm{0},\quad\bm{K}(\bm{U})=\bm{F}(\bm{U})-\bm{W}^{x}(\bm{U}),\quad\bm{L}(\bm{U})=\bm{G}(\bm{U})-\bm{W}^{y}(\bm{U}),$
(4.2)

where

$\bm{W}^{x}(\bm{U}):=\int\limits^{x}_{\widehat{x}}\left[B(\bm{U})\bm{U}_{\xi}(\xi,t)+\bm{S}^{x}(\bm{U}(\xi,t))\right]{\rm
d}\xi,\quad\bm{W}^{y}(\bm{U}):=\int\limits^{y}_{\widehat{y}}\left[C(\bm{U})\bm{U}_{\eta}(\eta,t)+\bm{S}^{y}(\bm{U}(\eta,t))\right]{\rm
d}\eta,$

with $\widehat{x}$ and $\widehat{y}$ being arbitrary numbers.

Applying the 2-D version of the semi-discrete LCD-CU scheme from §2 to the
quasi-conservative system (4.2) results in the following system of ODEs:

$\frac{{\rm d}\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j,k}}{{\rm
d}t}=-\frac{\bm{\mathcal{K}}^{\rm
LCD}_{{j+\frac{1}{2}},k}-\bm{\mathcal{K}}^{\rm
LCD}_{{j-\frac{1}{2}},k}}{\Delta x}-\frac{\bm{\mathcal{L}}^{\rm
LCD}_{j,{k+\frac{1}{2}}}-\bm{\mathcal{L}}^{\rm
LCD}_{j,{k-\frac{1}{2}}}}{\Delta y},$ (4.3)

where the numerical fluxes $\bm{\mathcal{K}}^{\rm LCD}_{{j+\frac{1}{2}},k}$
and $\bm{\mathcal{L}}^{\rm LCD}_{j,{k+\frac{1}{2}}}$ are given by

$\displaystyle\bm{\mathcal{K}}^{\rm LCD}_{{j+\frac{1}{2}},k}$
$\displaystyle=R_{{j+\frac{1}{2}},k}P^{\rm
LCD}_{{j+\frac{1}{2}},k}R^{-1}_{{j+\frac{1}{2}},k}\bm{K}^{-}_{{j+\frac{1}{2}},k}+R_{{j+\frac{1}{2}},k}M^{\rm
LCD}_{{j+\frac{1}{2}},k}R^{-1}_{{j+\frac{1}{2}},k}\bm{K}^{+}_{{j+\frac{1}{2}},k}$
(4.4) $\displaystyle+R_{{j+\frac{1}{2}},k}Q^{\rm
LCD}_{{j+\frac{1}{2}},k}R^{-1}_{{j+\frac{1}{2}},k}\big{(}\breve{\bm{U}}^{+}_{{j+\frac{1}{2}},k}-\breve{\bm{U}}^{-}_{{j+\frac{1}{2}},k}\big{)},$
$\displaystyle\bm{\mathcal{L}}^{\rm LCD}_{j,{k+\frac{1}{2}}}$
$\displaystyle=R_{j,{k+\frac{1}{2}}}P^{\rm
LCD}_{j,{k+\frac{1}{2}}}R^{-1}_{j,{k+\frac{1}{2}}}\bm{L}^{-}_{j,{k+\frac{1}{2}}}+R_{j,{k+\frac{1}{2}}}M^{\rm
LCD}_{j,{k+\frac{1}{2}}}R^{-1}_{j,{k+\frac{1}{2}}}\bm{L}^{+}_{j,{k+\frac{1}{2}}}$
$\displaystyle+R_{j,{k+\frac{1}{2}}}Q^{\rm
LCD}_{j,{k+\frac{1}{2}}}R^{-1}_{j,{k+\frac{1}{2}}}\big{(}\breve{\bm{U}}^{+}_{j,{k+\frac{1}{2}}}-\breve{\bm{U}}^{-}_{j,{k+\frac{1}{2}}}\big{)}.$

All of the terms in (4.4) are computed as in §3 separately in the $x$\- and
$y$-directions.

The global fluxes $\bm{K}^{\pm}_{{j+\frac{1}{2}},k}$ and
$\bm{L}^{\pm}_{j,{k+\frac{1}{2}}}$ are computed using the WB path-conservative
technique as in [2]. The LCD matrices $P^{\rm LCD}_{{j+\frac{1}{2}},k}$,
$M^{\rm LCD}_{{j+\frac{1}{2}},k}$, $Q^{\rm LCD}_{{j+\frac{1}{2}},k}$, $P^{\rm
LCD}_{j,{k+\frac{1}{2}}}$, $M^{\rm LCD}_{j,{k+\frac{1}{2}}}$, and $Q^{\rm
LCD}_{j,{k+\frac{1}{2}}}$ are computed as in [12], but based on the
eigensystems of $\widehat{\cal A}^{\,\bm{F}}_{{j+\frac{1}{2}},k}={\cal
A}^{\bm{F}}(\widehat{\bm{U}}_{{j+\frac{1}{2}},k}):=\frac{\partial\bm{F}}{\partial\bm{U}}(\widehat{\bm{U}}_{{j+\frac{1}{2}},k})-B(\widehat{\bm{U}}_{{j+\frac{1}{2}},k})$
and $\widehat{\cal A}^{\,\bm{G}}_{j,{k+\frac{1}{2}}}={\cal
A}^{\bm{G}}(\widehat{\bm{U}}_{j,{k+\frac{1}{2}}}):=\frac{\partial\bm{G}}{\partial\bm{U}}(\widehat{\bm{U}}_{j,{k+\frac{1}{2}}})-C(\widehat{\bm{U}}_{j,{k+\frac{1}{2}}})$.
As in §2, we do not employ local characteristic variables in the piecewise
linear reconstruction and perform the desingularization of $P^{\rm
LCD}_{{j+\frac{1}{2}},k}$, $M^{\rm LCD}_{{j+\frac{1}{2}},k}$, $Q^{\rm
LCD}_{{j+\frac{1}{2}},k}$, $P^{\rm LCD}_{j,{k+\frac{1}{2}}}$, $M^{\rm
LCD}_{j,{k+\frac{1}{2}}}$, and $Q^{\rm LCD}_{j,{k+\frac{1}{2}}}$ as in (2.8).

###### Remark 4.1

The introduced 2-D semi-discrete flux globalization based LCD-PCCU scheme is
WB in the sense that it is capable of exactly preserving discrete steady
states, for which $\bm{K}(\bm{U})$ is independent of $x$ and $\bm{L}(\bm{U})$
is independent of $y$. We omit the proof for the sake of brevity.

##  5 Numerical Examples

In this section, we apply the developed LCD-PCCU schemes to the $\gamma$-based
compressible multifluid and TRSW equations. We conduct several numerical
experiments and compare the performance of the LCD-PCCU schemes with their
PCCU counterparts obtained by applying the original CU numerical fluxes in
(3.1) and (4.3). For instance, in the 1-D case, such a PCCU scheme is obtained
by replacing (2.8) by

$\big{(}(P_{i})_{j+\frac{1}{2}},(M_{i})_{j+\frac{1}{2}},(Q_{i})_{j+\frac{1}{2}}\big{)}=\left\\{\begin{aligned}
&\frac{1}{a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}}\Big{(}a^{+}_{j+\frac{1}{2}},-a^{-}_{j+\frac{1}{2}},a^{+}_{j+\frac{1}{2}}a^{-}_{j+\frac{1}{2}}\Big{)}&&\mbox{if}~{}a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}>\varepsilon_{0},\\\\[1.29167pt]
&\Big{(}{\frac{1}{2}},{\frac{1}{2}},0\Big{)}&&\mbox{otherwise}.\end{aligned}\right.$

In all of the numerical examples, we have solved the ODE systems (3.1) and
(4.3) using the three-stage third-order strong stability preserving (SSP)
Runge-Kutta method (see, e.g., [23, 22]) and used the CFL number 0.45.

###  5.1 $\gamma$-Based Compressible Multifluid Equations

In this section, we show the applications of the developed LCD-PCCU schemes
introduced in §3 and §4 to the 1-D and 2-D $\gamma$-based compressible
multifluid systems [11, 39], which in the 2-D case can be written as (4.1)
with $\bm{U}=\big{(}\rho,\rho u,\rho v,E,\Gamma,\Pi\big{)}^{\top}$,
$\bm{F}(\bm{U})=\big{(}\rho u,\rho u^{2}+p,\rho uv,u(E+p),0,0\big{)}^{\top}$,
$\bm{G}(\bm{U})=\big{(}\rho v,\rho uv,\rho v^{2}+p,v(E+p),0,0)^{\top}$,
$\bm{S}^{x}(\bm{U})=\bm{S}^{y}(\bm{U})=\bm{0}$, $B(\bm{U})={\rm
diag}(0,0,0,0,u,u)$, and $C(\bm{U})={\rm diag}(0,0,0,0,v,v)$. Here, $\rho$ is
the density, $u$ and $v$ are the $x$\- and $y$-velocities, $p$ is the
pressure, and $E$ is the total energy. The system is closed through the
following (stiff) equation of state (EOS) for each of the fluid components:

$p=(\gamma-1)\left[E-\frac{\rho}{2}(u^{2}+v^{2})\right]-\gamma\pi_{\infty},$
(5.1)

where the parameters $\gamma$ and $\pi_{\infty}$ represent the specific heat
ratio and stiffness parameter, respectively (when $\pi_{\infty}\equiv 0$ for
all of the components, the system reduces to the ideal gas multicomponent
case). The variables $\Gamma$ and $\Pi$ are defined as
$\Gamma:=\frac{1}{\gamma-1}$ and $\Pi:=\frac{\gamma\pi_{\infty}}{\gamma-1}$.

###### Remark 5.1

We note that the studied $\gamma$-based compressible multifluid equations
contain no source terms and thus there is no need to numerically enforce the
WB property. We therefore do not need to introduce any equilibrium
variables—instead, we reconstruct the conservative quantities $\bm{U}$ and
replace $\breve{\bm{U}}^{\pm}_{j+\frac{1}{2}}$,
$\breve{\bm{U}}^{\pm}_{{j+\frac{1}{2}},k}$, and
$\breve{\bm{U}}^{\pm}_{j,{k+\frac{1}{2}}}$ in (3.2) and (4.4) by
$\bm{U}^{\pm}_{j+\frac{1}{2}}$, $\bm{U}^{\pm}_{{j+\frac{1}{2}},k}$, and
$\bm{U}^{\pm}_{j,{k+\frac{1}{2}}}$, respectively.

###### Remark 5.2

In Appendix A, we provide the expressions for the matrices $\widehat{\cal
A}_{j+\frac{1}{2}}$ and $\widehat{\cal A}^{\,\bm{F}}_{{j+\frac{1}{2}},k}$ and
the corresponding matrices $R_{j+\frac{1}{2}}$ and $R_{{j+\frac{1}{2}},k}$ for
the 1-D and 2-D $\gamma$-based compressible multifluid systems, respectively.
The matrices $\widehat{\cal A}^{\,\bm{G}}_{j,{k+\frac{1}{2}}}$ and
$R_{j,{k+\frac{1}{2}}}$ can be obtained in an analogous way.

####  5.1.1 One-Dimensional Examples

In this section, we consider two numerical examples for the 1-D version of the
$\gamma$-based compressible multifluid equations.

#### Example 1—“Shock-Bubble” Interaction

In the first example taken from [17], we consider the “shock-bubble”
interaction problem. The initial data are given by

$(\rho,u,p;\gamma,\pi_{\infty})(x,0)=\begin{cases}(13.1538,0;1,5/3,0),&|x|<0.25,\\\
(1.3333,-0.3535,1.5;1.4,0),&x>0.75,\\\
(1,0,1;1.4,0),&\mbox{otherwise,}\end{cases}$

which correspond to a left-moving shock, initially located at $x=0.75$, and a
resting “bubble” with a radius of 0.25, initially located at the origin.

We impose the solid wall boundary conditions at the both ends of the
computational domain $[-1,2]$ and compute the numerical solutions by both the
LCD-PCCU and PCCU schemes until the final time $t=3$ on a uniform mesh with
$\Delta x=1/100$. The obtained densities are presented in Figure 5.1 together
with the reference solution computed by the PCCU scheme on a much finer mesh
with $\Delta x=1/1000$. As one can see, the LCD-PCCU scheme achieves sharper
resolution than the PCCU one.

Figure 5.1: Example 1: Density $\rho$ computed by the PCCU and LCD-PCCU
schemes (left) and zoom at $x\in[-0.55,-0.4]$ (right).

#### Example 2—Water-Air “Shock-Bubble” Interaction

In the second 1-D example also taken from [17], we consider a gas-liquid
multifluid system, where the liquid component is modeled by the stiff EOS
(5.1) with $\pi_{\infty}\gg 1$. The initial conditions,

$(\rho,u,p;\gamma,\pi_{\infty})=\begin{cases}(0.05,0,1;1.4,0),&|x-6|<3,\\\
(1.325,-68.525,19153;4.4,6000),&x>11.4,\\\
(1,0,1;4.4,6000),&\mbox{otherwise},\end{cases}$

correspond to the left-moving shock, initially located at $x=11.4$, and a
resting air “bubble” with a radius 3, initially located at $x=6$.

We impose the free boundary conditions at the both ends of the computational
domain $[0,18]$ and compute the numerical solutions until the final time
$t=0.045$ on a uniform mesh with $\Delta x=1/10$ by both the LCD-PCCU and PCCU
schemes. The obtained results are presented in Figure 5.2 together with the
reference solution computed by the PCCU scheme on a much finer mesh with
$\Delta x=1/400$. As one can see, the LCD-PCCU solution achieves sharper
resolution than its PCCU counterpart, especially of the contact wave located
at about $x=3$.

Figure 5.2: Example 2: Density $\rho$ computed by the PCCU and LCD-PCCU
schemes (left) and zoom at $x\in[2.5,4.8]$ (right).

####  5.1.2 Two-Dimensional Examples

We now consider three 2-D examples, in which we plot the Schlieren images of
the magnitude of the density gradient field, $|\nabla\rho|$. To this end, we
use the shading function

$\exp\bigg{(}-\frac{80|\nabla\rho|}{\max(|\nabla\rho|)}\bigg{)},$

where the numerical derivatives of the density are computed using standard
central differencing.

#### Example 3—Shock-Helium Bubble Interaction

In the first 2-D example taken from [37, 15], a shock wave in the air hits the
light resting bubble which contains helium. We take the following initial
conditions:

$(\rho,u,v,p;\gamma,p_{\infty})=\begin{cases}(4/29,0,0,1;5/3,0)&\mbox{in
region A},\\\ (1,0,0,1;1.4,0)&\mbox{in region B},\\\
(4/3,-0.3535,0,1.5;1.4,0)&\mbox{in region C},\end{cases}$

where regions A, B, and C are outlined in Figure 5.3, and the computational
domain is $[-3,1]\times[-0.5,0.5]$. We impose the solid wall boundary
conditions on the top and bottom and the free boundary conditions on the left
and right edges of the computational domain.

Figure 5.3: Initial setting for the 2-D numerical examples.

We compute the numerical solutions until the final time $t=3$ on a uniform
mesh with $\Delta x=\Delta y=1/500$. In Figure 5.4, we present different
stages of the shock-bubble interaction computed by the PCCU and LCD-PCCU
schemes. One can clearly see that the bubble interface develops very complex
structures after the bubble is hit by the shock, especially at large times. We
note that the obtained results are in a good agreement with the experimental
findings presented in [24] and the computational results presented in [37, 15,
13]. Notably, at an early time $t=0.5$, the resolution of the bubble interface
significantly improves with the implementation of the LCD-PCCU scheme. As time
progresses, instabilities arise within the interface, which are smeared by a
more dissipative PCCU scheme. The differences in resolution become more
pronounced at later times $t=2$ and $2.5$, and especially at $t=3$, indicating
that the LCD-PCCU scheme clearly outperforms the PCCU one.

Figure 5.4: Example 3: Shock-helium bubble interaction by the PCCU (left
column) and LCD-PCCU (right column) schemes at different times.

#### Example 4—Shock-R22 Bubble Interaction

In the second 2-D example also taken from [37, 15], a shock wave in the air
hits the heavy resting bubble which contains R22 gas. The initial conditions
are

$(\rho,u,v,p;\gamma,\pi_{\infty})=\begin{cases}(3.1538,0,0,1;1.249,0)&\mbox{in
region A},\\\ (1,0,0,1;1.4,0)&\mbox{in region B},\\\
(4/3,-0.3535,0,1.5;1.4,0)&\mbox{in region C},\end{cases}$

and the rest of the settings are the same as in Example 3.

We compute the numerical solutions until the final time $t=3$ on a uniform
mesh with $\Delta x=\Delta y=1/500$ and present different stages of the shock-
bubble interaction computed by the PCCU and LCD-PCCU schemes in Figure 5.5.
One can clearly see that upon encountering the bubble, the shock wave
undergoes both refraction and reflection. Under the force of the shock, the
bubble undergoes compression and motion. These findings are in a good
agreement with the computational results reported in [37, 15]. Figure 5.5
demonstrates that the LCD-PCCU scheme captures material interfaces more
accurately (compared to the PCCU scheme) and also enhances the resolution of
fine details of the computed solution. Once again, this indicates that the
LCD-PCCU scheme is less diffusive than its PCCU counterpart.

Figure 5.5: Example 4: Shock-R22 bubble interaction by the PCCU (left column)
and LCD-PCCU (right column) schemes at different times.

#### Example 5—A Cylindrical Explosion Problem

In the last 2-D example taken from [17] (see also [41]), we consider the case
where a cylindrical explosive source is located between an air-water interface
and an impermeable wall. The initial conditions are given by

$(\rho,u,v,p;\gamma,\pi_{\infty})=\begin{cases}(1.27,0,0,8290,2,0),&\mbox{$(x-5)^{2}+(y-2)^{2}<1$},\\\
(0.02,0,0,1;1.4,0),&\mbox{$y>4$},\\\
(1,0,0,1;7.15,3309),&\mbox{otherwise},\end{cases}$

with the solid wall boundary conditions imposed at the bottom, and the free
boundary conditions set on the other sides of the computational domain
$[0,10]\times[0,6]$.

We compute the numerical solutions until the final time $t=0.02$ on a uniform
mesh with $\Delta x=\Delta y=1/80$ by the studied PCCU and LCD-PCCU schemes.
In Figure 5.6, we present time snapshots of the obtained results, which
qualitatively look similar to the numerical results reported in [17, 41]. As
one can see, the LCD-PCCU scheme captures both the material interfaces and
many of the developed wave structures substantially sharper than the PCCU
scheme.

Figure 5.6: Example 5: Time evolution of the solution computed by the PCCU
(left column) and LCD-PCCU (right column) schemes.

One can also notice a substantial difference between the PCCU and LCD-PCCU
solutions at the final time. In order to clarify this point, we refine the
mesh and compute the numerical solution on a finer mesh with $\Delta x=\Delta
y=1/160$ by the PCCU scheme and present the obtained result in Figure 5.7,
where one can clearly see that the PCCU solution computed on a finer mesh is
close to the LCD-PCCU solution computed with $\Delta x=\Delta y=1/80$. Once
again, this indicates that the LCD-PCCU solution is more accurate than the
PCCU one.

Figure 5.7: Example 5: A finer mesh ($\Delta x=\Delta y=1/160$) PCCU solution
at $t=0.02$.

###  5.2 Thermal Rotating Shallow Water Equations

In this section, we apply the developed LCD-PCCU schemes to the 1-D and 2-D
TRSW equations [27], which in the 2-D case can be written as the system (4.1)
with

$\begin{aligned}
&\bm{U}=(h,q,p,hb)^{\top},\quad\bm{F}(\bm{U})=\Big{(}hu,hu^{2}+\frac{b}{2}h^{2},huv,hub\Big{)}^{\top},\quad\bm{G}(\bm{U})=\Big{(}hv,huv,hv^{2}+\frac{b}{2}h^{2},hvb\Big{)}^{\top},\\\
&B(\bm{U})=C(\bm{U})\equiv\bm{0},\quad\bm{S}^{x}(\bm{U})=(0,fhv-
hbZ_{x},0,0)^{\top},\quad\bm{S}^{y}(\bm{U})=(0,0,-fhu-
hbZ_{y},0)^{\top},\end{aligned}$

(5.2)

where $h$ denotes the thickness of the fluid layer, $u$ and $v$ are the zonal
and meridional velocities, $q:=hu$ and $p:=hv$ are the corresponding
discharges, $b$ is the layer-averaged buoyancy variable, and $Z(x,y)$ stands
for a time-independent bottom topography. In the oceanic context,
$b=g\rho/\rho_{0}$, where $g$ is the acceleration due to gravity, $\rho$ and
$\rho_{0}$ are the variable and constant parts of the water density,
respectively. In the atmospheric context, the densities $\rho$ and $\rho_{0}$
should be replaced with the potential temperatures $\theta$ and $\theta_{0}$.
Finally, $f(y)=f_{0}+\beta y$ is the Coriolis parameter, where $f_{0}$ and
$\beta$ are positive constants.

###### Remark 5.3

For the TRSW equations, we use the equilibrium variables introduced in [2] and
compute the point values ($\bm{U}^{\pm}_{k+\frac{1}{2}}$ and
$\breve{\bm{U}}^{\pm}_{k+\frac{1}{2}}$ in the 1-D case and
$\bm{U}^{\pm}_{{j+\frac{1}{2}},k}$, $\bm{U}^{\pm}_{j,{k+\frac{1}{2}}}$,
$\breve{\bm{U}}^{\pm}_{{j+\frac{1}{2}},k}$, and
$\breve{\bm{U}}^{\pm}_{j,{k+\frac{1}{2}}}$ in the 2-D case) as well as the
global fluxes ($\bm{L}^{\pm}_{k+\frac{1}{2}}$ in the 1-D case and
$\bm{K}^{\pm}_{{j+\frac{1}{2}},k}$ and $\bm{L}^{\pm}_{j,{k+\frac{1}{2}}}$ in
the 2-D case) precisely as it was done in [2]. We have also used the same
numerical diffusion switch function, which helps to enforce the WB property,
as in [2].

###### Remark 5.4

In Appendix B, we provide the expressions for the matrices $\widehat{\cal
A}_{k+\frac{1}{2}}$ and $\widehat{\cal A}^{\,\bm{F}}_{{j+\frac{1}{2}},k}$ and
the corresponding matrices $R_{k+\frac{1}{2}}$ and $R_{{j+\frac{1}{2}},k}$ for
the 1-D and 2-D TRSW equations, respectively. The matrices $\widehat{\cal
A}^{\,\bm{G}}_{j,{k+\frac{1}{2}}}$ and $R_{j,{k+\frac{1}{2}}}$ can be obtained
in an analogous way.

####  5.2.1 One-Dimensional Examples

In this section, we consider two numerical examples for the so-called 1-D Ripa
system, which is the 1-D version of (4.1), (5.2), but with $f(y)\equiv 0$:

$\begin{cases}h_{t}+(hv)_{y}=0,\\\\[1.29167pt]
p_{t}+\Big{(}hv^{2}+\dfrac{b}{2}h^{2}\Big{)}_{y}=-hbZ_{y},\\\\[2.15277pt]
(hb)_{t}+(hvb)_{y}=0.\end{cases}$

#### Example 6—Constant Pressure Equilibrium

In the first example taken from [2], we consider the following initial data:

$(h,v,b)(y,0)=\begin{cases}(2,0,1)&\mbox{if}~{}y<0,\\\
(1,0,4)&\mbox{otherwise},\end{cases}$ (5.3)

which are prescribed in the computational domain $[-5,5]$ together with the
flat bottom topography $Z(y)\equiv 0$ subject to the free boundary conditions.

We first compute the numerical solution by the PCCU and LCD-PCCU schemes until
the final time $t=10$ on the uniform mesh with $\Delta x=1/20$ and plot the
obtained fluid thickness and buoyancy variables in Figure 5.8, where the
reference solution is computed by the PCCU scheme on a much finer mesh with
$\Delta x=1/400$. One can clearly see that both the PCCU and LCD-PCCU schemes
can exactly preserve the steady state.

Figure 5.8: Example 6: $h$ and $b$ computed by the PCCU and LCD-PCCU schemes.

We then test the ability of the proposed schemes to capture small
perturbations of this steady state. To this end, we first denote by $h_{\rm
eq}(y)=h(y,0)$ given by (5.3), and then modify the initial condition of $h$ in
the following manner:

$h(y,0)=h_{\rm eq}(y)+\begin{cases}0.1&\mbox{if}~{}y\in[-1.8,-1.7],\\\
0&\mbox{otherwise}.\end{cases}$

We compute the solution until the final time $t=1.6$ on a uniform mesh with
$\Delta x=1/20$ and plot the obtained fluid thickness at different times
$t=1.2$ and 1.6 in Figure 5.9 together with the reference solution computed by
the LCD-PCCU scheme on a much finer mesh with $\Delta x=1/400$. One can
clearly see that the LCD-PCCU scheme achieves sharper resolution than the PCCU
scheme. Moreover, the LCD-PCCU solution is substantially less oscillatory than
the PCCU one.

Figure 5.9: Example 6: $h$ computed by the PCCU and LCD-PCCU schemes at
$t=1.2$ (top row) and 1.6 (bottom row) and zoom at $x\in[-1.5,1]$ (right
column).

#### Example 7—Dam Break Over a Nonflat Bottom

In the second example taken from [16], we consider the dam break problem
subject to the following initial data:

$(h,v,b)(y,0)=\begin{cases}(5-Z(y),0,1)&\mbox{if}~{}y<0,\\\
(2-Z(y),0,5)&\mbox{otherwise},\end{cases}$

and the nonflat bottom topography

$Z(y)=\begin{cases}\cos\big{(}10\pi(x+0.3)\big{)}+1&\mbox{if}~{}-0.4\leq
y\leq-0.2,\\\
{\frac{1}{2}}\big{[}\cos\big{(}10\pi(x-0.3)\big{)}+1\big{]}&0.2\leq y\leq
0.4,\\\ 0&\mbox{otherwise}.\end{cases}$

The initial conditions are prescribed in the computational domain $[-1,1]$
subject to the free boundary conditions.

We compute the numerical solution by the PCCU and LCD-PCCU schemes until the
final time $t=0.2$ on a uniform mesh with $\Delta x=1/100$ and plot the
obtained water surface $h+Z$ in Figure 5.10, where the reference solution is
computed by the PCCU scheme on a much finer mesh with $\Delta x=1/2000$. One
can clearly see that the LCD-PCCU scheme produces higher-resolution results
than the PCCU one.

Figure 5.10: Example 7: $h+Z$ computed by the PCCU and LCD-PCCU schemes at
$t=0.1$ (top row) and 0.2 (bottom row) and zoom at $x\in[-0.05,0.1]$ (right
column).

####  5.2.2 Two-Dimensional Examples

In this section, we consider three numerical examples for the 2-D TRSW
equations.

#### Example 8—Small Perturbations of a 2-D Steady State

In this example taken from [16], we take $f(y)\equiv 0$ and consider the
following initial conditions:

$(h,u,v,b)(x,y,0)=\begin{cases}\big{(}3.1-Z(x,y),0,0,\frac{4}{3}\big{)}&\mbox{if}~{}0.01<x^{2}+y^{2}<0.09,\\\
\big{(}3-Z(x,y),0,0,\frac{4}{3}\big{)}&\mbox{if}~{}0.09<x^{2}+y^{2}<0.25~{}{\rm
or}~{}x^{2}+y^{2}<0.01,\\\ (2-Z(x,y),0,0,3)&\mbox{otherwise},\end{cases}$

with the nonflat bottom topography consisting of two Gaussian shaped humps:

$Z(x,y)=\begin{cases}0.5\,e^{-100\left[(x+0.5)^{2}+(y+0.5)^{2}\right]}&\mbox{if}~{}x<0,\\\
0.6\,e^{-100\left[(x-0.5)^{2}+(y-0.5)^{2}\right]}&\mbox{otherwise}.\end{cases}$

The initial data are prescribed in the computational domain
$[-1,1]\times[-1,1]$ subject to the free boundary conditions.

We compute the numerical solution by the PCCU and LCD-PCCU schemes until the
final time $t=0.12$ on a uniform mesh with $\Delta x=\Delta y=1/50$. The
obtained results are presented in Figure 5.11, where one can clearly see that
the LCD-PCCU scheme outperforms the PCCU one, since the $h$\- and
$b$-components computed by the PCCU scheme are smeared; see Figure 5.11. We
also plot the 1-D slices of $h$ and $b$ along the diagonal $y=x$; see Figure
5.12, where the reference solution is computed by the PCCU scheme on a much
finer mesh with $\Delta x=\Delta y=1/400$. As one can see, the LCD-PCCU scheme
achieves higher resolution than the PCCU one.

Figure 5.11: Example 8: $h$ and $b$ computed by the PCCU (left column) and
LCD-PCCU (right column) schemes.

Figure 5.12: Example 8: Diagonal slices of $h$ (left) and $b$ (right) computed
by the PCCU and LCD-PCCU schemes; zoom at $[-0.75,0.75]$.

#### Example 9—Radial Dam Break Over the Flat Bottom

In this example taken from [16], we compare the PCCU and LCD-PCCU schemes on a
circular dam-break problem with a flat bottom topography ($Z(x,y)\equiv 0$)
and no Coriolis forces $(f(y)\equiv 0)$. The initial conditions,

$(h,u,v,b)(x,y,0)=\begin{cases}(1.5,0,0,1)&\mbox{if}~{}x^{2}+y^{2}<0.25,\\\
(1.2,0,0,1.5)&\mbox{otherwise},\end{cases}$

are prescribed in the computational domain $[-1,1]\times[-1,1]$ subject to the
free boundary conditions.

We compute the numerical solution by the PCCU and LCD-PCCU schemes until the
final time $t=0.15$ on a uniform mesh with $\Delta x=\Delta y=1/50$. The
computed results are presented in Figure 5.13, where one can clearly see that
the LCD-PCCU scheme outperforms the PCCU scheme. We also plot the 1-D slices
of $h$ and $b$ along the diagonal $y=x$; see Figure 5.14, where the reference
solution is computed by the PCCU scheme on a much finer mesh with $\Delta
x=\Delta y=1/400$. As one can see, the LCD-PCCU scheme outperforms its PCCU
counterpart.

Figure 5.13: Example 9: $h$ and $b$ computed by the PCCU (left column) and
LCD-PCCU (right column) schemes.

Figure 5.14: Example 9: Diagonal slices of $h$ (left) and $b$ (right) computed
by the PCCU and LCD-PCCU schemes; zoom at $[-0.6,0.6]$.

#### Example 10—Relaxation of Localized Anomalies on the Equatorial
$\beta$-Plane

In the last example taken from [27], we consider the relaxation of localized
anomalies on the $\beta$-plane, where the Coriolis parameter is $f(y)=y$, with
the initial data

$h(x,y,0)\equiv 1,\quad q(x,y,0)=p(x,y,0)\equiv 0,\quad
b(x,y,0)=1+{\frac{1}{2}}\,e^{-\left(\frac{1}{50}x^{2}+\frac{1}{2}y^{2}\right)},$

prescribed in the computational domain $[-40,80]\times[-10,10]$ subject to the
free boundary conditions.

We compute the numerical solution by the PCCU and LCD-PCCU schemes until the
final time $t=120$ on a uniform mesh with $\Delta x=\Delta y=2/15$ and plot
the obtained numerical results in Figure 5.15, where one can see that our
results are consistent with those reported in [27], and at the same time, the
LCD-PCCU scheme outperforms the PCCU one. In order to better see the
difference between the computed solutions, we also plot the 1-D slices of $h$
along $y=0$ in Figure 5.16, where the reference solution is computed by the
PCCU scheme on a finer mesh with $\Delta x=\Delta y=1/25$. One can see that
the LCD-PCCU scheme achieves significantly better resolution than the PCCU
one.

Figure 5.15: Example 10: Time snapshots of $h$ computed by the PCCU (left
column) and LCD-PCCU (right column) schemes.

Figure 5.16: Example 10: Slices of $h$ and $b$ along the line $y=0$ computed
by the PCCU and LCD-PCCU schemes using $\Delta x=\Delta y=2/15$ at different
times.

In this example, we also check the efficiency of the proposed LCD-PCCU scheme
by performing a more thorough comparison between the studied schemes. To this
end, we measure the CPU time consumed during the above calculations by the
LCD-PCCU scheme and refine the mesh used by the PCCU scheme to the level that
exactly the same CPU time is consumed to compute both of the numerical
solutions: This way, we take into account additional computational cost of the
LCD-PCCU scheme. The corresponding grids are $\Delta x=\Delta y=2/15$ for the
LCD-PCCU scheme, and $\Delta x=\Delta y=1/9$ for the PCCU scheme. The obtained
numerical results are plotted in Figure 5.17, where one can clearly see that
the LCD-PCCU solution is still more accurate than the PCCU one.

Figure 5.17: Example 10: Slices of $h$ and $b$ along the line $y=0$ computed
by the PCCU scheme using $\Delta x=\Delta y=1/9$ and LCD-PCCU scheme using
$\Delta x=\Delta y=2/15$ at different times.

### Acknowledgments

The work of Michael Herty was funded by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation)–SPP 2183: Eigenschaftsgeregelte
Umformprozesse with the Project(s) HE5386/19-2,19-3 Entwicklung eines
flexiblen isothermen Reckschmiedeprozesses für die eigenschaftsgeregelte
Herstellung von Turbinenschaufeln aus Hochtemperaturwerkstoffen (424334423)
and by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation)–SPP 2410 Hyperbolic Balance Laws in Fluid Mechanics: Complexity,
Scales, Randomness (CoScaRa) within the Project(s) HE5386/26-1 (Numerische
Verfahren für gekoppelte Mehrskalenprobleme,525842915) and (Zufällige
kompressible Euler Gleichungen: Numerik und ihre Analysis, 525853336)
HE5386/27-1. The work of A. Kurganov was supported in part by NSFC grant
12171226 and by the fund of the Guangdong Provincial Key Laboratory of
Computational Science and Material Design (No. 2019B030301001).

##  Appendix A LCD Matrices for the $\gamma$-Based Compressible Multifluid
Equations

###  A.1 One-Dimensional Case

In order to compute the LCD matrices for the 1-D $\gamma$-based compressible
multifluid system, we first compute the matrix ${\cal
A}=\frac{\partial\bm{F}}{\partial\bm{U}}-B$:

${\cal A}(\bm{U})=\begin{pmatrix}0&1&0&0&0\\\
\dfrac{\gamma-3}{2}u^{2}&(3-\gamma)u&\gamma-1&(1-\gamma)p&1-\gamma\\\\[3.44444pt]
\dfrac{uc^{2}}{1-\gamma}+\Big{(}\dfrac{\gamma}{2}-1\Big{)}u^{3}&\dfrac{c^{2}}{\gamma-1}+\Big{(}\dfrac{3}{2}-\gamma\Big{)}u^{2}&\gamma
u&(1-\gamma)pu&(1-\gamma)u\\\ 0&0&0&u&0\\\ 0&0&0&0&u\end{pmatrix},$ (A.1)

where $c=\sqrt{\frac{\gamma(p+\pi_{\infty})}{\rho}}$ is the speed of sound,
and then introduce $\widehat{\cal A}_{j+\frac{1}{2}}={\cal
A}(\widehat{\bm{U}}_{j+\frac{1}{2}})$, which is given by (A.1) with $\gamma$,
$u$, $p$, and $c$ replaced with

$\hat{\gamma}=\frac{\gamma_{j}+\gamma_{j+1}}{2},\quad\hat{u}=\frac{u_{j}+u_{j+1}}{2},\quad\hat{p}=\frac{p_{j}+p_{j+1}}{2},\quad{\rm
and}\quad\hat{c}=\sqrt{\frac{\hat{\gamma}(\hat{p}+\hat{\pi}_{\infty})}{\hat{\rho}}},$

respectively. Here,

$\displaystyle\gamma_{j}=\frac{1}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}_{j}}+1,\quad(\pi_{\infty})_{j}=\frac{\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern
0.0pt}}}_{j}}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}_{j}+1},\quad\hat{\rho}=\frac{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}_{j}+\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern
0.0pt}}}_{j+1}}{2},\quad u_{j}=\frac{(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})_{j}}{\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}_{j}},$
$\displaystyle p_{j}=(\gamma_{j}-1)\Big{[}\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern
0.0pt}}}_{j}-{\frac{1}{2}}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern
0.0pt}}}_{j}u_{j}^{2}\Big{]}-\gamma_{j}(\pi_{\infty})_{j},\quad\hat{\pi}_{\infty}=\frac{(\pi_{\infty})_{j}+(\pi_{\infty})_{j+1}}{2}.$

Notice that all of the $\hat{(\cdot)}$ quantities have to have a subscript
index, that is, $\hat{(\cdot)}=\hat{(\cdot)}_{j+\frac{1}{2}}$, but we omit it
for the sake of brevity for all of the quantities except for $\widehat{\cal
A}_{j+\frac{1}{2}}$. We then compute the matrix $R_{j+\frac{1}{2}}$ composed
of the right eigenvectors of $\widehat{\cal A}_{j+\frac{1}{2}}$ and obtain

$\displaystyle R_{j+\frac{1}{2}}=$ $\displaystyle\begin{pmatrix}1&1&0&0&1\\\
\hat{u}-\hat{c}&\hat{u}&0&0&\hat{u}+\hat{c}\\\\[1.29167pt]
\dfrac{\hat{c}^{2}}{\hat{\gamma}-1}+\dfrac{1}{2}\hat{u}^{2}-\hat{u}\hat{c}&\dfrac{\hat{u}^{2}}{2}&\hat{p}&0&\dfrac{\hat{c}^{2}}{\hat{\gamma}-1}+\dfrac{1}{2}\hat{u}^{2}+\hat{u}\hat{c}\\\\[2.15277pt]
0&0&1&-1&0\\\ 0&0&0&\hat{p}&0\end{pmatrix},$ $\displaystyle
R^{-1}_{j+\frac{1}{2}}=\dfrac{1}{2\hat{c}^{2}}$
$\displaystyle\begin{pmatrix}\dfrac{(\hat{\gamma}-1)\hat{u}^{2}}{2}+\hat{u}\hat{c}&-\hat{c}-(\hat{\gamma}-1)\hat{u}&\hat{\gamma}-1&(1-\hat{\gamma})\hat{p}&1-\hat{\gamma}\\\\[2.15277pt]
2\hat{c}^{2}-(\hat{\gamma}-1)\hat{u}^{2}&2(\hat{\gamma}-1)\hat{u}&2(1-\hat{\gamma})&2(\hat{\gamma}-1)\hat{p}&2(\hat{\gamma}-1)\\\\[1.29167pt]
0&0&0&2\hat{c}^{2}&\dfrac{2\hat{c}^{2}}{\hat{p}}\\\\[6.45831pt]
0&0&0&0&-\dfrac{2\hat{c}^{2}}{\hat{p}}\\\
\dfrac{(\hat{\gamma}-1)\hat{u}^{2}}{2}-\hat{u}\hat{c}&\hat{c}-(\hat{\gamma}-1)\hat{u}&\hat{\gamma}-1&(1-\hat{\gamma})\hat{p}&1-\hat{\gamma}\end{pmatrix}.$

###  A.2 Two-Dimensional Case

We now compute the matrix $\widehat{\cal A}^{\,\bm{F}}_{{j+\frac{1}{2}},k}$
and the corresponding matrices $R_{{j+\frac{1}{2}},k}$ and
$R^{-1}_{{j+\frac{1}{2}},k}$ for the 2-D $\gamma$-based compressible
multifluid system.

We first compute ${\cal A}^{\bm{F}}=\frac{\partial\bm{F}}{\partial\bm{U}}-B$:

${\cal A}^{\bm{F}}(\bm{U})=\begin{pmatrix}0&1&0&0&0&0\\\
\dfrac{\gamma-3}{2}u^{2}+\dfrac{\gamma-1}{2}v^{2}&(3-\gamma)u&(1-\gamma)v&\gamma-1&(1-\gamma)p&1-\gamma\\\\[2.15277pt]
-uv&v&u&0&0&0\\\
\dfrac{uc^{2}}{1-\gamma}+\Big{(}\dfrac{\gamma}{2}-1\Big{)}u(u^{2}+v^{2})&\dfrac{c^{2}}{\gamma-1}+\Big{(}\dfrac{3}{2}-\gamma\Big{)}u^{2}+\dfrac{1}{2}v^{2}&(1-\gamma)uv&\gamma
u&(1-\gamma)pu&(1-\gamma)u\\\ 0&0&0&0&u&0\\\ 0&0&0&0&0&u\end{pmatrix},$

(A.2)

and then introduce the matrix $\widehat{\cal
A}^{\,\bm{F}}_{{j+\frac{1}{2}},k}={\cal
A}(\widehat{\bm{U}}_{{j+\frac{1}{2}},k})$, which is given by (A.2) with
$\gamma$, $u$, $v$, $p$, and $c$ replaced with

$\hat{\gamma}=\frac{\gamma_{j,k}+\gamma_{j+1,k}}{2},\quad\hat{u}=\frac{u_{j,k}+u_{j+1,k}}{2},\quad\hat{v}=\frac{v_{j,k}+v_{j+1,k}}{2},\quad\hat{p}=\frac{p_{j,k}+p_{j+1,k}}{2},\quad\hat{c}=\sqrt{\frac{\hat{\gamma}(\hat{p}+\hat{\pi}_{\infty})}{\hat{\rho}}},$

respectively. Here,

$\displaystyle\gamma_{j,k}=\frac{1}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}_{j,k}}+1,\quad(\pi_{\infty})_{j,k}=\frac{\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern
0.0pt}}}_{j,k}}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}_{j,k}+1},\quad\hat{\rho}=\frac{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}_{j,k}+\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern
0.0pt}}}_{j+1,k}}{2},\quad u_{j,k}=\frac{(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern
0.0pt}}})_{j,k}}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}_{j,k}},\quad
v_{j,k}=\frac{(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern
0.0pt}}})_{j,k}}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}_{j,k}},$ $\displaystyle
p_{j,k}=(\gamma_{j,k}-1)\Big{[}\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern
0.0pt}}}_{j,k}-{\frac{1}{2}}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern
0.0pt}}}_{j,k}u_{j,k}^{2}-{\frac{1}{2}}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern
0.0pt}}}_{j,k}v_{j,k}^{2}\Big{]}-\gamma_{j,k}(\pi_{\infty})_{j,k},\quad\hat{\pi}_{\infty}=\frac{(\pi_{\infty})_{j,k}+(\pi_{\infty})_{j+1,k}}{2}.$

Notice that all of the $\hat{(\cdot)}$ quantities have to have a subscript
index, that is, $\hat{(\cdot)}=\hat{(\cdot)}_{{j+\frac{1}{2}},k}$, but we omit
it for the sake of brevity for all of the quantities except for $\widehat{\cal
A}^{\,\bm{F}}_{{j+\frac{1}{2}},k}$. We then compute the matrix
$R_{{j+\frac{1}{2}},k}$ composed of the right eigenvectors of $\widehat{\cal
A}^{\,\bm{F}}_{{j+\frac{1}{2}},k}$ and obtain

$\displaystyle R_{{j+\frac{1}{2}},k}=$
$\displaystyle\begin{pmatrix}1&1&0&0&0&1\\\
\hat{u}-\hat{c}&\hat{u}&0&0&0&\hat{u}+\hat{c}\\\
\hat{v}&\hat{v}&1&0&0&\hat{v}\\\
\dfrac{\hat{c}^{2}}{\hat{\gamma}-1}+\dfrac{1}{2}\hat{u}^{2}+\dfrac{1}{2}\hat{v}^{2}-\hat{u}\hat{c}&\dfrac{\hat{u}^{2}+\hat{v}^{2}}{2}&\hat{v}&\hat{p}&0&\dfrac{\hat{c}^{2}}{\hat{\gamma}-1}+\dfrac{1}{2}\hat{u}^{2}+\dfrac{1}{2}\hat{v}^{2}+\hat{u}\hat{c}\\\
0&0&0&1&1&0\\\ 0&0&0&0&-\hat{p}&0\end{pmatrix}$ $\displaystyle
R^{-1}_{{j+\frac{1}{2}},k}=\dfrac{1}{2\hat{c}^{2}}$
$\displaystyle\begin{pmatrix}\dfrac{(\hat{\gamma}-1)(\hat{u}^{2}+\hat{v}^{2})}{2}+\hat{u}\hat{c}&-\hat{c}-(\hat{\gamma}-1)\hat{u}&(1-\hat{\gamma})\hat{v}&\hat{\gamma}-1&(1-\hat{\gamma})\hat{p}&1-\hat{\gamma}\\\\[2.15277pt]
2\hat{c}^{2}-(\hat{\gamma}-1)(\hat{u}^{2}+\hat{v}^{2})&2(\hat{\gamma}-1)\hat{u}&2(\hat{\gamma}-1)\hat{v}&2(1-\hat{\gamma})&2(\hat{\gamma}-1)\hat{p}&2(\hat{\gamma}-1)\\\
-2\hat{c}^{2}\hat{v}&0&2\hat{c}^{2}&0&0&0\\\\[1.29167pt]
0&0&0&0&2\hat{c}^{2}&\dfrac{2\hat{c}^{2}}{\hat{p}}\\\
\dfrac{(\hat{\gamma}-1)(\hat{u}^{2}+\hat{v}^{2})}{2}-\hat{u}\hat{c}&\hat{c}-(\hat{\gamma}-1)hatu&(1-\hat{\gamma})\hat{v}&\hat{\gamma}-1&(1-\hat{\gamma})\hat{p}&1-\hat{\gamma}\end{pmatrix}.$

##  Appendix B LCD Matrices for the TRSW Equations

###  B.1 One-dimensional Case

In order to compute the LCD matrices for the 1-D TRSW equations, we first
compute the matrix ${\cal A}=\frac{\partial\bm{G}}{\partial\bm{U}}-C$:

${\cal A}(\bm{U})=\begin{pmatrix}0&0&1&0\\\ -uv&v&u&0\\\\[1.29167pt]
\dfrac{bh}{2}-v^{2}&0&2v&\dfrac{h}{2}\\\\[1.29167pt] -bv&0&b&v\end{pmatrix},$
(B.1)

and then introduce the matrices $\widehat{\cal A}_{k+\frac{1}{2}}={\cal
A}(\widehat{\bm{U}}_{k+\frac{1}{2}})$, which are given by (B.1) with $h$, $u$,
$v$, and $b$ replaced with

$\hat{h}=\frac{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$h$\kern 0.0pt}}}_{k}+\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h$\kern
0.0pt}}}_{k+1}}{2},\quad\hat{u}=\frac{u_{k}+u_{k+1}}{2},\quad\hat{v}=\frac{v_{k}+v_{k+1}}{2},\quad\hat{b}=\frac{b_{k}+b_{k+1}}{2},$

respectively. Here, $u_{k}=\nicefrac{{\hbox{\vbox{\hrule height=0.5pt\kern
1.20552pt\hbox{\kern-0.35002pt$q$\kern 0.0pt}}}_{k}}}{{\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.20552pt\hbox{\kern-0.35002pt$h$\kern 0.0pt}}}_{k}}}$,
$v_{k}=\nicefrac{{\hbox{\vbox{\hrule height=0.5pt\kern
1.20552pt\hbox{\kern-0.35002pt$p$\kern 0.0pt}}}_{k}}}{{\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.20552pt\hbox{\kern-0.35002pt$h$\kern 0.0pt}}}_{k}}}$, and
$b_{k}=\nicefrac{{\hbox{\vbox{\hrule height=0.5pt\kern
1.20552pt\hbox{\kern-0.35002pt$(hb)$\kern 0.0pt}}}_{k}}}{{\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.20552pt\hbox{\kern-0.35002pt$h$\kern 0.0pt}}}_{k}}}$.

Notice that all of the $\hat{(\cdot)}$ quantities have to have a subscript
index, that is, $\hat{(\cdot)}=\hat{(\cdot)}_{k+\frac{1}{2}}$, but we omit it
for the sake of brevity for all of the quantities except for $\widehat{\cal
A}_{k+\frac{1}{2}}$. We then compute the matrix $R_{k+\frac{1}{2}}$ composed
of the right eigenvectors of $\widehat{\cal A}_{k+\frac{1}{2}}$ and obtain

$R_{k+\frac{1}{2}}=\frac{1}{\hat{b}}\begin{pmatrix}1&-1&0&1\\\
\hat{u}&0&\hat{b}&\hat{u}\\\
\hat{v}-\sqrt{\hat{b}\hat{h}}&-\hat{v}&0&\hat{v}+\sqrt{\hat{b}\hat{h}}\\\
\hat{b}&\hat{b}&0&\hat{b}\end{pmatrix},\quad
R^{-1}_{k+\frac{1}{2}}=\frac{1}{4}\begin{pmatrix}\hat{b}+2\hat{v}\sqrt{\dfrac{\hat{b}}{\hat{h}}}&0&-2\sqrt{\dfrac{\hat{b}}{\hat{h}}}&1\\\\[7.74998pt]
-2\hat{b}&0&0&2\\\ -2\hat{u}&4&0&-2\dfrac{\hat{u}}{\hat{b}}\\\\[4.30554pt]
\hat{b}-2\hat{v}\sqrt{\dfrac{\hat{b}}{\hat{h}}}&0&2\sqrt{\dfrac{\hat{b}}{\hat{h}}}&1\end{pmatrix}.$

###  B.2 Two-Dimensional Case

We now compute the matrix $\widehat{\cal A}^{\,\bm{F}}_{{j+\frac{1}{2}},k}$
and the corresponding matrices $R_{{j+\frac{1}{2}},k}$ and
$R^{-1}_{{j+\frac{1}{2}},k}$ for the 2-D TRSW equations.

We first compute ${\cal
A}^{\,\bm{F}}=\frac{\partial\bm{F}}{\partial\bm{U}}-B$:

${\cal A}^{\,\bm{F}}(\bm{U})=\begin{pmatrix}0&1&0&0\\\
\dfrac{bh}{2}-u^{2}&2u&0&\dfrac{h}{2}\\\\[3.44444pt] -uv&v&u&0\\\
-bu&b&0&u\end{pmatrix},$ (B.2)

and then introduce the matrices $\widehat{\cal
A}^{\,\bm{F}}_{{j+\frac{1}{2}},k}={\cal
A}(\widehat{\bm{U}}_{{j+\frac{1}{2}},k})$, which are given by (B.2) with $h$,
$u$, $v$, and $b$ replaced with

$\hat{h}=\frac{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$h$\kern 0.0pt}}}_{j,k}+\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h$\kern
0.0pt}}}_{j+1,k}}{2},\quad\hat{u}=\frac{u_{j,k}+u_{j+1,k}}{2},\quad\hat{v}=\frac{v_{j,k}+v_{j+1,k}}{2},\quad\hat{b}=\frac{b_{j,k}+b_{j+1,k}}{2},$

respectively. Here $u_{k}=\nicefrac{{\hbox{\vbox{\hrule height=0.5pt\kern
1.20552pt\hbox{\kern-0.35002pt$q$\kern 0.0pt}}}_{j,k}}}{{\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.20552pt\hbox{\kern-0.35002pt$h$\kern 0.0pt}}}_{j,k}}}$,
$v_{k}=\nicefrac{{\hbox{\vbox{\hrule height=0.5pt\kern
1.20552pt\hbox{\kern-0.35002pt$p$\kern 0.0pt}}}_{j,k}}}{{\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.20552pt\hbox{\kern-0.35002pt$h$\kern 0.0pt}}}_{j,k}}}$,
and $b_{k}=\nicefrac{{\hbox{\vbox{\hrule height=0.5pt\kern
1.20552pt\hbox{\kern-0.35002pt$(hb)$\kern
0.0pt}}}_{j,k}}}{{\,\hbox{\vbox{\hrule height=0.5pt\kern
1.20552pt\hbox{\kern-0.35002pt$h$\kern 0.0pt}}}_{j,k}}}$.

Notice that all of the $\hat{(\cdot)}$ quantities have to have a subscript
index, that is, $\hat{(\cdot)}=\hat{(\cdot)}_{{j+\frac{1}{2}},k}$, but we omit
it for the sake of brevity for all of the quantities except for $\widehat{\cal
A}^{\,\bm{F}}_{{j+\frac{1}{2}},k}$. We then compute the matrix
$R_{j+\frac{1}{2}}$ composed of the right eigenvectors of $\widehat{\cal
A}^{\,\bm{F}}_{{j+\frac{1}{2}},k}$ and obtain

$R_{{j+\frac{1}{2}},k}=\frac{1}{\hat{b}}\begin{pmatrix}1&-1&0&1\\\
\hat{u}-\sqrt{\hat{b}\hat{h}}&-\hat{u}&0&\hat{u}+\sqrt{\hat{b}\hat{h}}\\\
\hat{v}&0&\hat{b}&\hat{v}\\\ \hat{b}&\hat{b}&0&\hat{b}\end{pmatrix},\quad
R^{-1}_{{j+\frac{1}{2}},k}=\frac{1}{4}\begin{pmatrix}\hat{b}+2\hat{u}\sqrt{\dfrac{\hat{b}}{\hat{h}}}&0&-2\sqrt{\dfrac{\hat{b}}{\hat{h}}}&1\\\\[7.74998pt]
-2\hat{b}&0&0&2\\\ -2\hat{v}&0&4&-2\dfrac{\hat{v}}{\hat{b}}\\\\[4.30554pt]
\hat{b}-2\hat{u}\sqrt{\dfrac{\hat{b}}{\hat{h}}}&2\sqrt{\dfrac{\hat{b}}{\hat{h}}}&0&1\end{pmatrix}.$

## References

  * [1] S. Busto, M. Dumbser, S. Gavrilyuk, and K. Ivanova, On thermodynamically compatible finite volume methods and path-conservative ADER discontinuous Galerkin schemes for turbulent shallow water flows, J. Sci. Comput., 88 (2021).  Paper No. 28, 45 pp.
  * [2] Y. Cao, A. Kurganov, and Y. Liu, Flux globalization based well-balanced path-conservative central-upwind scheme for the thermal rotating shallow water equations, Commun. Comput. Phys., 34 (2023), pp. 993–1042.
  * [3] Y. Cao, A. Kurganov, Y. Liu, and R. Xin, Flux globalization based well-balanced path-conservative central-upwind schemes for shallow water models, J. Sci. Comput., 92 (2022).  Paper No. 69, 31 pp.
  * [4] Y. Cao, A. Kurganov, Y. Liu, and V. Zeitlin, Flux globalization based well-balanced path-conservative central-upwind scheme for two-layer thermal rotating shallow water equations, J. Comput. Phys., 474 (2023).  Paper No. 111790, 29 pp.
  * [5] V. Caselles, R. Donat, and G. Haro, Flux-gradient and source-term balancing for certain high resolution shock-capturing schemes, Comput. & Fluids, 38 (2009), pp. 16–36.
  * [6] M. J. Castro, T. Morales de Luna, and C. Parés, Well-balanced schemes and path-conservative numerical methods, in Handbook of numerical methods for hyperbolic problems, vol. 18 of Handb. Numer. Anal., Elsevier/North-Holland, Amsterdam, 2017, pp. 131–175.
  * [7] M. J. Castro and C. Parés, Well-balanced high-order finite volume methods for systems of balance laws, J. Sci. Comput., 82 (2020).  Paper No. 48, 48 pp.
  * [8] M. J. Castro Díaz, A. Kurganov, and T. Morales de Luna, Path-conservative central-upwind schemes for nonconservative hyperbolic systems, ESAIM Math. Model. Numer. Anal., 53 (2019), pp. 959–985.
  * [9] C. Chalons, Path-conservative in-cell discontinuous reconstruction schemes for non conservative hyperbolic systems, Commun. Math. Sci., 18 (2020), pp. 1–30.
  * [10] Y. Chen, A. Kurganov, and M. Na, A Flux Globalization Based Well-Balanced Path-Conservative Central-Upwind Scheme for the Shallow Water Flows in Channels, ESAIM Math. Model. Numer. Anal., 57 (2023), pp. 1087–1110.
  * [11] J. Cheng, F. Zhang, and T. G. Liu, A discontinuous Galerkin method for the simulation of compressible gas-gas and gas-water two-medium flows, J. Comput. Phys., 403 (2020).  Paper No. 109059, 29 pp.
  * [12] A. Chertock, S. Chu, M. Herty, A. Kurganov, and M. Lukáčová-Medviďová, Local characteristic decomposition based central-upwind scheme, J. Comput. Phys., 473 (2023).  Paper No. 111718, 24 pp.
  * [13] A. Chertock, S. Chu, and A. Kurganov, Hybrid multifluid algorithms based on the path-conservative central-upwind scheme, J. Sci. Comput., 89 (2021).  Paper No. 48, 24 pp.
  * [14] A. Chertock, M. Herty, and Ş. N. Özcan, Well-balanced central-upwind schemes for $2\times 2$ systems of balance laws, in Theory, Numerics and Applications of Hyperbolic Problems. I, C. Klingenberg and M. Westdickenberg, eds., vol. 236 of Springer Proc. Math. Stat., Springer, Cham, 2018, pp. 345–361.
  * [15] A. Chertock, S. Karni, and A. Kurganov, Interface tracking method for compressible multifluids, M2AN Math. Model. Numer. Anal., 42 (2008), pp. 991–1019.
  * [16] A. Chertock, A. Kurganov, and Y. Liu, Central-upwind schemes for the system of shallow water equations with horizontal temperature gradients, Numer. Math., 127 (2014), pp. 595–639.
  * [17] S. Chu, A. Kurganov, and R. Xin, Low-dissipation central-upwind schemes for compressible multifluids, Submitted.  Preprint available at https://sites.google.com/view/alexander-kurganov/publications.
  * [18] S. Chu, A. Kurganov, and R. Xin, A well-balanced fifth-order A-WENO scheme based on flux globalization, Submitted.  Preprint available at https://sites.google.com/view/alexander-kurganov/publications.
  * [19] G. Dal Maso, P. G. Lefloch, and F. Murat, Definition and weak stability of nonconservative products, J. Math. Pures Appl., 74 (1995), pp. 483–548.
  * [20] D. Donat and A. Martinez-Gavara, Hybrid second order schemes for scalar balance laws, J. Sci. Comput., 48 (2011), pp. 52–69.
  * [21] L. Gascón and J. M. Corderán, Construction of second-order TVD schemes for nonhomogeneous hyperbolic conservation laws, J. Comput. Phys., 172 (2001), pp. 261–297.
  * [22] S. Gottlieb, D. Ketcheson, and C.-W. Shu, Strong stability preserving Runge-Kutta and multistep time discretizations, World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2011.
  * [23] S. Gottlieb, C.-W. Shu, and E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM Rev., 43 (2001), pp. 89–112.
  * [24] J.-F. Haas and B. Sturtevant, Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities, J. Fluid Mech., 181 (1987), pp. 41–76.
  * [25] A. Kurganov and C.-T. Lin, On the reduction of numerical dissipation in central-upwind schemes, Commun. Comput. Phys., 2 (2007), pp. 141–163.
  * [26] A. Kurganov, Y. Liu, and R. Xin, Well-balanced path-conservative central-upwind schemes based on flux globalization, J. Comput. Phys., 474 (2023).  Paper No. 111773, 32 pp.
  * [27] A. Kurganov, Y. Liu, and V. Zeitlin, Thermal versus isothermal rotating shallow water equations: comparison of dynamical processes by simulations with a novel well-balanced central-upwind scheme, Geophys. Astrophys. Fluid Dyn., 115 (2021), pp. 125–154.
  * [28] A. Kurganov, S. Noelle, and G. Petrova, Semidiscrete central-upwind schemes for hyperbolic conservation laws and Hamilton-Jacobi equations, SIAM J. Sci. Comput., 23 (2001), pp. 707–740.
  * [29] A. Kurganov and E. Tadmor, New high-resolution semi-discrete central schemes for Hamilton-Jacobi equations, J. Comput. Phys., 160 (2000), pp. 720–742.
  * [30] Philippe G. LeFloch, Hyperbolic systems of conservation laws, Lectures in Mathematics ETH Zürich, Birkhäuser Verlag, Basel, 2002.  The theory of classical and nonclassical shock waves.
  * [31] P. G. Lefloch, Graph solutions of nonlinear hyperbolic systems, J. Hyperbolic Differ. Equ., 1 (2004), pp. 643–689.
  * [32] K.-A. Lie and S. Noelle, On the artificial compression method for second-order nonoscillatory central difference schemes for systems of conservation laws, SIAM J. Sci. Comput., 24 (2003), pp. 1157–1174.
  * [33] A. Martinez-Gavara and R. Donat, A hybrid second order scheme for shallow water flows, J. Sci. Comput., 48 (2011), pp. 241–257.
  * [34] H. Nessyahu and E. Tadmor, Nonoscillatory central differencing for hyperbolic conservation laws, J. Comput. Phys., 87 (1990), pp. 408–463.
  * [35] C. Parés, Path-conservative numerical methods for nonconservative hyperbolic systems, vol. 24 of Quad. Mat., Dept. Math., Seconda Univ. Napoli, Caserta, 2009.
  * [36] E. Pimentel-García, M. J. Castro, C. Chalons, T. Morales de Luna, and C. Parés, In-cell discontinuous reconstruction path-conservative methods for non conservative hyperbolic systems—second-order extension, J. Comput. Phys., 459 (2022).  Paper No. 111152, 35 pp.
  * [37] J. J. Quirk and S. Karni, On the dynamics of a shock-bubble interaction., J. Fluid Mech., 318 (1996), pp. 129–163.
  * [38] K. A. Schneider, J. M. Gallardo, D. S. Balsara, B. Nkonga, and C. Parés, Multidimensional approximate Riemann solvers for hyperbolic nonconservative systems. Applications to shallow water systems, J. Comput. Phys., 444 (2021).  Paper No. 110547, 49 pp.
  * [39] K.-M. Shyue, An efficient shock-capturing algorithm for compressible multicomponent problems, J. Comput. Phys., 142 (1998), pp. 208–242.
  * [40] P. K. Sweby, High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM J. Numer. Anal., 21 (1984), pp. 995–1011.
  * [41] L. Xu and T. G. Liu, Explicit interface treatments for compressible gas-liquid simulations, Comput. & Fluids, 153 (2017), pp. 34–48.