Papers by Piotr Smolarkiewicz
Scientific Final Report: COLLABORATIVE RESEARCH: CONTINUOUS DYNAMIC GRID ADAPTATION IN A GLOBAL ATMOSPHERIC MODEL: APPLICATION AND REFINEMENT
This project had goals of advancing the performance capabilities of the numerical general circula... more This project had goals of advancing the performance capabilities of the numerical general circulation model EULAG and using it to produce a fully operational atmospheric global climate model (AGCM) that can employ either static or dynamic grid stretching for targeted phenomena. The resulting AGCM combined EULAG's advanced dynamics core with the 'physics' of the NCAR Community Atmospheric Model (CAM). Effort discussed below shows how we improved model performance and tested both EULAG and the coupled CAM-EULAG in several ways to demonstrate the grid stretching and ability to simulate very well a wide range of scales, that is, multi-scale capability. We leveraged our effort through interaction with an international EULAG community that has collectively developed new features and applications of EULAG, which we exploited for our own work summarized here. Overall, the work contributed to over 40 peer-reviewed publications and over 70 conference/workshop/seminar presentations, many of them invited.
Key numerical issues for future development of the ECMWF models
... numerical issues for the future development of the ECMWF model By Mike Cullen, Deborah Salmon... more ... numerical issues for the future development of the ECMWF model By Mike Cullen, Deborah Salmond and ... The step which enforces the conti-nuity equation is still elliptic in nature, even though in ... Even if the problem lo be solved is nonlinear, only a linear problem can be solved ...
Gravity waves (GWs) are ubiquitous internal waves in which the restoring force is buoyancy. In th... more Gravity waves (GWs) are ubiquitous internal waves in which the restoring force is buoyancy. In the Earth's atmosphere, GWs are generated whenever a parcel of air is disturbed from its equilibrium position. Disturbing sources include topographic features (e.g., mountains), convective instabilities (e.g., thunderstorms), frontal motions, and oscillations of jet streams. When the amplitude of a GW approaches its inverse vertical wavenumber, the wave "breaks" and generates a localized region of intense turbulence. Strong GWs near the surface in the lee of mountains are responsible for spectacular down-slope windstorms, such as the Boulder chinook, the Santa Ana winds, and the Yugoslavian bora. Such storms often generate hurricane-force winds and are a serious threat to personal safety, buildings, and surface and air transportation. At higher altitudes, GWs can be of sufficient magnitude to alter the mean state of the atmosphere.
Journal of the Atmospheric Sciences, Sep 1, 1996
An anelastic approximation is used with a time-variable coordinate transformation to formulate a ... more An anelastic approximation is used with a time-variable coordinate transformation to formulate a two-dimensional numerical model that describes the evolution of gravity waves. The model is solved using a semi-Lagrangian method with monotone (nonoscillatory) interpolation of all advected fields. The time-variable transformation is used to generate disturbances at the lower boundary that approximate the effect of a traveling line of thunderstorms (a squall line) or of flow over a broad topographic obstacle. The vertical propagation and breaking of the gravity wave field (under conditions typical of summer solstice) is illustrated for each of these cases. It is shown that the wave field at high altitudes is dominated by a single horizontal wavelength, which is not always related simply to the horizontal dimension of the source. The morphology of wave breaking depends on the horizontal wavelength; for sufficiently short waves, breaking involves roughly one half of the wavelength. In common with other studies, it is found that the breaking waves undergo `self-acceleration,' such that the Zonal-mean intrinsic frequency remains approximately constant in spite of large changes in the background wind. It is also shown that many of the features obtained in the calculations can be understood in terms of linear wave theory. In particular, linear theory provides insights into the wavelength of the waves that break at high altitudes, the onset and evolution of breaking, the horizontal extent of the breaking region and its position relative to the forcing, and the minimum and maximum altitudes where breaking occurs. Wave breaking ceases at the altitude where the background dissipation rate (which in our model is a proxy for molecular diffusion) becomes greater than the rate of dissipation due to wave breaking. This altitude, in effect, the model turbopause, is shown to depend on a relatively small number of parameters that characterize the waves and the background state.
A viscoelastic fluid model for brain injuries
International Journal For Numerical Methods in Fluids, Sep 9, 2002
SUMMARY Due to its elasticity, the human brain material can support shear (equivoluminal) waves. ... more SUMMARY Due to its elasticity, the human brain material can support shear (equivoluminal) waves. Earlier attempts to explain certain brain injuries via arguments of the classical theory of viscoelasticity exploited the Voigt model—a linear system of dierential equations where the motion of the brain tissue depends merely on the balance between viscous and elastic forces. Although Voigt model solutions illustrate the role of the viscoelastic mechanics in brain injuries, they have limited use for modelling realistic cases which, for example, evince strongly localized displacements of the brain tissue. We have extended the Voigt model to a non-linear viscoelasticuid model, thereby dispensing with simplifying assumptions of vanishing advective transport. The resulting non-Newtonianuid model admits non-linear phenom- ena such as steepening of the wave fronts as well as wave overturning and their subsequent turbulent breaking. The posed equations are solved numerically, and the solution procedure are validated against small-perturbation linear theory and closed-form Voigt-model solutions available in the literature. Our non-linear numerical results suggest existence of a 'brain turbulence' phenomenon. They are in qualita- tive agreement with the results of medical research, especially, with regard to the diuse axonal injuries which are observed to occur in a highly localized manner near the border between the gray and the white matter. Copyright ? 2002 John Wiley & Sons, Ltd.
The Effects of Forcing Modes on the Simulations of Tidally Generated Internal Bores and Solitons
Numerical experiments were performed on the generation of internal waves by barotropic tidal curr... more Numerical experiments were performed on the generation of internal waves by barotropic tidal currents flowing over strongly varying topography, such as sills in straits. The studies concentrated on two aspects of the simulations: (a) the amplitude of the generated waves as functions of tidal velocities, and (b) amplitudes obtained with different modes of forcing prescriptions, such as boundary vs. volume forcing, or analytic vs. model-generated currents. The tidal velocities were obtained from three tidal components, the M2, K1 and O1 tides, and the model currents were computed with a hydrostatic model. In the boundary-forcing approach, the current variations are prescribed as normal inflow at the respective boundaries. In the volume-forcing approach, the currents are prescribed at all interior points, with the model currents already adapted to the interior topography, but the analytic sinusoidal currents had to undergo a flux-conservation modification. Simulations over both idealized topographies, and the sills of the Luzon Strait and South China Sea were carried out. Propagation comparisons were also made with those of the "shape-preserving", analytically constructed type solitons (solutions of the KDV equation).. For tidal currents in the 10-30 cm/sec range, the generated perturbation densities typically ranged peak-to-peak from 1.6 to 3.5 kg/m3, and the vertical velocities from 25 to 70 cm/sec. The analytic solitons maintained their amplitudes during propagation better than the tidally-generated signals.
Applied Mathematics and Computation, Jun 1, 2001
A portable Fortran code that solves the general nonseparable three-dimensional linear elliptic pa... more A portable Fortran code that solves the general nonseparable three-dimensional linear elliptic partial dierential equation (PDE) with cross-derivative terms on a rectangular region is described. Boundary conditions can be any combination of periodic, speci®ed, or mixed derivative. A multigrid scheme, modi®ed to handle complexities introduced by cross-derivative terms and nonnormal derivative components in boundary conditions, is utilized. Successful application is illustrated with an analytically prescribed abstract example and simulation of an idealized geophysical¯uid ow. Ó
We present a systematic study of numerical aspects of advection/diffusion in thermal convection a... more We present a systematic study of numerical aspects of advection/diffusion in thermal convection as realized in large-scale high-resolution numerical models. With rapid progress in computer technology the NWP models with 10 km horizontal resolution (or better) will become standard in the near future; cf. recent calculations by the Frontier Research Center group in Japan (Satoh et al, J. Comput. Phys. ) conducted at 3.5 km resolution. While such resolutions are impressive by the standards of NWP, they are still overly coarse (by two orders of magnitude at least) to represent convection up to the standards of cloudresolving models. In effect, NWP is entering a new regime, where traditional convection parameterization is already obsolete, but large-eddy-simulation is still beyond the reach. In this regime the convective response to advection/diffusion realization is extremely sensitive to filtering embedded in numerical model -e.g., via subgrid-scale models and/or numerical approximations employed -and can take variety of forms. This in turn reflects on simulated weather and climate due to their dependence on cloud field structure via precipitation and radiation. Here, we report on our investigation of numerical effects that influence structure of simulated convection at its roots; i.e., in planetary boundary layer. In particular, we focus on effective viscosity and diffusivity representative of contemporary numerical models.
Scientific Final Report: COLLABORATIVE RESEARCH: CONTINUOUS DYNAMIC GRID ADAPTATION IN A GLOBAL ATMOSPHERIC MODEL: APPLICATION AND REFINEMENT
This project had goals of advancing the performance capabilities of the numerical general circula... more This project had goals of advancing the performance capabilities of the numerical general circulation model EULAG and using it to produce a fully operational atmospheric global climate model (AGCM) that can employ either static or dynamic grid stretching for targeted phenomena. The resulting AGCM combined EULAG's advanced dynamics core with the 'physics' of the NCAR Community Atmospheric Model (CAM). Effort discussed below shows how we improved model performance and tested both EULAG and the coupled CAM-EULAG in several ways to demonstrate the grid stretching and ability to simulate very well a wide range of scales, that is, multi-scale capability. We leveraged our effort through interaction with an international EULAG community that has collectively developed new features and applications of EULAG, which we exploited for our own work summarized here. Overall, the work contributed to over 40 peer-reviewed publications and over 70 conference/workshop/seminar presentations, many of them invited.
Coupling the dynamics of boundary layers and evolutionary dunes
Physical Review E Statistical Nonlinear and Soft Matter Physics, May 1, 2009
A theoretical formulation and corresponding numerical solutions are presented for fluid flow and ... more A theoretical formulation and corresponding numerical solutions are presented for fluid flow and sediment transport past evolutionary sand dunes. Time-dependent curvilinear coordinates are employed to fully couple flow aloft with the developing landform. The differential conservation law that defines shape of the lower boundary depends on details of local surface stress, thereby favoring the large eddy simulation of the boundary layer. To shrink the gap between the time scales characteristic of planetary boundary layer flows O(103)s and sand dune evolution O(106)s , a hypothetical “severe-wind scenario” is adopted with the saltation flux amplified up to 3 orders of magnitude. While the results are largely insensitive to the rescaling, the efficacy of computations is greatly improved. The flux-form partial differential equation for the interface profile—via saltation and sand avalanches—is formulated as an advection-diffusion equation, to facilitate discrete integrations. Numerical experiments verify the adopted theoretical framework by reproducing scaling results reported in the literature. The versatility of the approach is illustrated with evolution of a sandhole—an example of application likely never addressed in the literature, yet realizable in nature.
Large-eddy simulations of convective boundary layers using nonoscillatory differencing
Physica D Nonlinear Phenomena, Sep 1, 1999
We explore the ability of nonoscillatory advection schemes to represent the effects of the unreso... more We explore the ability of nonoscillatory advection schemes to represent the effects of the unresolved scales of motion in numerical simulation of turbulent flows. We demonstrate that a nonoscillatory fluid solver can accurately reproduce the dynamics of an atmospheric convective boundary layer. When an explicit turbulence model is implemented, the solver does not add any significant numerical diffusion. Of greater interest, when no explicit turbulence is employed the solver itself appears to include an effective subgrid scale model. Other researchers have reported similar success, simulating turbulent flows in a variety of regimes while using only nonoscillatory advection schemes to model subgrid effects. At this point there is no theory justifying this success, but we offer some speculations.
Untwisting magnetic fields in the solar corona
ABSTRACT The solar corona is the tenuous atmosphere of the Sun characterized by a temperature of ... more ABSTRACT The solar corona is the tenuous atmosphere of the Sun characterized by a temperature of the order of million degrees Kelvin, an ambient magnetic field of 10 to 15 Gauss and a very high magnetic Reynolds number because of which it qualifies as a near-ideal magnetofluid system. It is well known that for such a system, the magnetic flux across every fluid surface remains effectively constant to a good approximation. Under this so called ``frozen-in'' condition then, it is possible to partition this magnetofluid into contiguous magnetic subvolumes each entrapping its own subsystem of magnetic flux. Thin magnetic flux tubes are an elementary example of such magnetic subvolumes evolving in time with no exchange of fluid among them. The internal twists and interweaving of these flux tubes, collectively referred as the magnetic topology, remains conserved under the frozen-in condition. Because of the dynamical evolution of the magnetofluid, two such subvolumes can come into direct contact with each other by expelling a third interstitial subvolume. In this process, the magnetic field may become discontinuous across the surface of contact by forming a current sheet there. Because of the small spatial scales generated by steepening of magnetic field gradient, the otherwise negligible resistivity becomes dominant and allows for reconnection of field lines which converts magnetic energy into heat. This phenomenon of spontaneous current sheet formation and its subsequent resistive decay is believed to be a possible mechanism for heating the solar corona to its million degree Kelvin temperature. In this work the dynamics of spontaneous current sheet formation is explored through numerical simulations and the results are presented.
Large scale numerical experiments on solar convection zone
37Th Cospar Scientific Assembly, 2008
Large scale dynamics under the influence of developing magnetic fields inside the solar convectio... more Large scale dynamics under the influence of developing magnetic fields inside the solar convection zone is investigated using the MHD version of EULAG, a global numerical model designed to work in the spirit of implicit large eddy simulations, using a higher-order upwind type advection scheme, in the framework of the anelastic approximation for MHD equations. Differential rotation profiles, turbulence characteristics and magnetic field evolution are analyzed under various parameters settings in an attempt to gain a better understanding of different factors contributions to the solar activity characteristics.
Cyclic dynamo action in numerial simulation of solar convection
ABSTRACT I will review recent advances in the production of solar-like cyclic variations of a dyn... more ABSTRACT I will review recent advances in the production of solar-like cyclic variations of a dynamo-generated large-scale magnetic field arising in global MHD simulations of solar convection. After describing a few representative simulations, I will focus on the turbulent nature of the driving electromotive force. Cyclic variations of the large-scale magnetic component appears to hinge on a fine balance between induction by the large-scale flows, and a turbulent electromotive force, the latter surprisingly akin to simple expectations from mean-field electrodynamics. I will also discuss the nature and mode of operatin of the mechanisms limiting the growth of the large-scale magnetic field in these simulations.
A 3D nonhydrostatic, Navier-Stokes solver has been employed to simulate gravity wave induced turb... more A 3D nonhydrostatic, Navier-Stokes solver has been employed to simulate gravity wave induced turbulence at mesopause altitudes. This paper extends our earlier 2D study reported in the literature to three spatial dimensions while maintaining fine resolution required to capture essential physics of the wave breaking. The calculations were performed on the 512 processor Cray T3E machine at the National Energy Research Scientific Computing Center (NERSC) in Berkeley. The physical results of this study clearly demonstrate advantages of highly parallel technologies. We briefly outline the physical outcome of the study, as well as compare the relative model performance across several machines using both MPI and Shmem communication software.
Ambient Hydrographic Effects on the Generation and Propagation of Tidally-Driven Internal Bores and Solitons
Agu Spring Meeting Abstracts, May 1, 2005
The effects of different initial hydrographic representations on the numerical simulations of int... more The effects of different initial hydrographic representations on the numerical simulations of internal bore and soliton generation and propagation were studied. The initial set of experiments were carried out in the Luzon Strait and the eastern part of the South China Sea using nonhydrostatic models. Internal wave generation was forced by the barotropic tide passing over the sills between islands and south of Taiwan in the Strait. Horizontally homogeneous initial density profiles were derived, via analytical approximations, from three sources: (1) the 2001 Levitus monthly, 1/4 deg climatology; (2) a set of 694 CTD observations collected between April and July of 2001, and (3) numerical simulations for the period of the observations. To obtain high resolution 3-D simulations, only subregions of the Luzon Strait/South China Sea area were considered, with focus on regions in the 118E-122E and 18.5N-22.5N geographic domain. Model resolutions varied from 50m to 1 km in the horizontal and 25m to 100m in the vertical. The choice of temperature and salinity profiles, and the resulting density profile, affected the numerical convergence of the iterations in the model. In addition, the amplitude of the generated internal bores and their propagation were also affected.
A Monge-Amp�re enhancement for semi-Lagrangian methods
Comput Fluids, 2011
Demanding the compatibility of semi-Lagrangian trajectory schemes with the fundamental Euler expa... more Demanding the compatibility of semi-Lagrangian trajectory schemes with the fundamental Euler expansion formula leads to the Monge-Ampère (MA) nonlinear second-order partial differential equation. Given standard estimates of the departure points of flow trajectories, solving the associated MA problem provides a corrected solution satisfying a discrete Lagrangian form of the mass continuity equation to round-off error. The impact of the MA enhancement is discussed in two diverse limits of fluid dynamics applications: passive tracer advection in a steady cellular flow and in fully developed turbulence. Improvements of the overall accuracy of simulations depend on the problem and can be substantial.
Intrigued by the regularity of convective structures observed in simulations of mesoscale flow pa... more Intrigued by the regularity of convective structures observed in simulations of mesoscale flow past realistic topography, we take a deeper look into computational aspects of a classical problem of the flow over a heated plane. We found that the numerical solutions are sensitive to viscosity, either incorporated a priori or effectively realized in computational models. In particular, anisotropic viscosity can lead to regular convective structures that mimic naturally occurring Rayleigh-Bénard (RB) cells but that are spurious for the problem at hand. We have extended the classical linear theory to anisotropic viscosity at moderately supercritical Rayleigh numbers, realized effectively in under-resolved convection simulations. It follows that anisotropic viscosity modifies the range of unstable RB modes, such that for an effective viscosity much larger in the horizontal than in the vertical unphysically broad RB cells may be observed. The latter is relevant to "cloud resolving" global models with relatively fine (for numerical weather prediction) horizontal resolution ^x ~O(103) m. At such a resolution the simulated convection is still under-resolved and strongly influenced by numerical filtering. To better assess the impact of an effective model viscosity on the structure of convective fields, we have conducted a large series of simulations of thermal convection, with various degrees of idealization, using the computational model EULAG. We performed an extensive convergence study and documented differences between the well resolved (viz. realistic) cellular convection and spurious structures. Comparing various means of enhancing the effective viscosity in the horizontal, we demonstrated that details of filtering are inessential. The common denominator of the scale selection is consistent with the linear theory. On the practical side we found that some numerical approaches may be preferable when the resolution is inadequate to capture the realism of convective fields. While control of effective viscosity is certainly the key to the quality results, resorting to non-dissipative numerics is not a cure. We found that implicit large-eddy simulation (ILES) approach based on non-oscillatory forward-in-time numerics minimizes numerical viscosity and its anisotropy, and produces results superior compared to more standard LES models.
Numerical Simulations of Closed Head Injuries
Spontaneous Current Sheet Formation and Break-Up of Magnetic Flux Surfaces
The dynamics of spontaneous current sheet formation is demonstrated in a viscous, perfectly condu... more The dynamics of spontaneous current sheet formation is demonstrated in a viscous, perfectly conducting, incompressible magnetofluid through numerical simulations. The magnetic field is represented in terms of evolving flux surfaces which are the possible sites of current sheet formation. The computation follows global magnetic flux surfaces of simple initial geometry as they deform into more complex forms creating current sheets in the process. Ultimately, the flux surfaces break their initial topology, as the spatial scale of surface folds decreases below the model resolution. This breaking is used to identify the sites of the current sheets formation.
Uploads
Papers by Piotr Smolarkiewicz