 Good morning, everyone. It's my great pleasure to be here. I've been working in the field of tectonic modeling and surface dynamics is not a total stranger in tectonic modeling field, but still I'm pretty much new to serious surface dynamics research. So I hope my title doesn't give you a false impression that I'm going to give a grand vision of unifying those two fields. Instead, I'll be a lot less ambitious and I'll take this opportunity to introduce to you a tool that I developed for tectonic modeling, which can potentially be useful for surface dynamics research. So here's the outline of my talk. First, I give an overview of the code snack, and then I'll give you some examples of tectonic models with surface processes included. Assuming that you are all familiar with erosion in organic environment, I chose a couple of examples from extensional environments, and at the end I'll give you an quick preview of ongoing project that is trying to couple one surface modeling code with the tectonic modeling code. And I'll discuss how snack can fit into that project. More often than not, I have to worry about losing audience because of too much math, but this CSDMS meeting must be one of the exceptions, and so I put up the equation slide first without being apologetic. I don't know, maybe I should apologize that there's too little math or too low-level math. Anyways, so snack is an acronym for St. Germain Analysis of Continua. St. Germain is the name of software framework on which snack is built. So name itself doesn't have much meaning, but LSD sounds friendly. The governing equation of snack is the momentum balance equation given in the form of the principle of virtual power. It is the use of this integrated form that rendered the snack comparable to finite element method. In snack, the spatial domain is discretized into hexahedron elements, and each hexahedron is subdivided into two different sets of five tetrahedron. So the basic unit of computation or the element type is the linear tetrahedron. But the mesh is organized in this hierarchical and redundant way to ensure a symmetric response. The final form of discretized momentum equation is nothing more than the Newton's second law of motion. So the time rate of change of momentum at a given global node L is the sum of surface force and body force, which are also the sum of all the contributing tetrahedron around that global node L. So once you get the net force acting on a node on a node, then acceleration is time-integrated to get velocity and update grid point. The one advantage of explicit formulation is that there is no need for nonlinear iteration. Because the time integration is explicit, all the updated quantities, velocity and grid point location, is depending only on the previous time steps quantities. And velocity and node locations are updated through leapfrog method. And it's obvious that description of motion is Lagrangian, in which grid points move with material. It's nice to have dynamic formulation because in the nonlinear problems we don't have to worry about convergence issue. However, to get a static or quasi-static solution, we need to damp out the inertia. That means we need to include a certain dissipative process. And SNAC uses so-called local damping, such that before we update the velocity, the actual net force on the node is damped by a certain damping force. And the damping force is determined as fixed proportion of the net force itself, but it's acting on the negative, the opposite direction to the current velocity. So if you look at this very simple 1D damped oscillator, it might sound ridiculous to go through this dynamic relaxation scheme to get a solution. Because you simply divide this gravitational Mg with spring constant K. But in multi-degrees of freedom problem, the spring constant K becomes very large sparse matrix, stiff matrix, and inverting it is not trivial at all. So this dynamic relaxation scheme is particularly useful for nonlinear problems. Another issue with dynamic formulation is that the time integration is not always stable, and the time step size is governed by so-called Quran condition. So if we have a kilometer resolution, and if we consider the typical P wave velocity in the crust, then time scale is easily becomes less than second. And that is prohibitive to cover millions or tens of million years of model time. One way to go around that limitation is to adjust the P wave velocity. So if we convert the usual definition of P wave velocity from stiffness matrix, sorry, stiffness divided by density, but if we volume integrate them, we can have corresponding stiffness divided by a certain mass. And we adjust this value of mass, we call it scaled mass, then we can have a lot smaller P wave velocity than physical value, then we can have sufficiently large time step. So the actual gravitational mass is not touched, so the problem itself doesn't change, only the dynamic relaxation step is affected. The constitutive model employed in SNAC is elasto-viscoplastic. And actually it's a linear combination of elastic, viscous, and plastic units. This model is very flexible. For instance, if we have a very high viscosity, then the material, whole material becomes essentially plastic material, for example, more Coulomb. And if yield stress is said to be extremely high, then whole material becomes essentially Maxwell viscoelastic. So if we have a temperature dependent viscosity, then burrito ductile transition will be naturally determined by a given temperature field. And if we have yield stress that is a function of a certain internal variable, then we can have strain hardening or softening, which is required for a strain localization phenomena. SNAC is an MPI parallel code. So here I'm showing some data from performance tests. So on the right, the plot shows work clock time for three different problem size as a function of number of cores. And on the right, it's a corresponding speed-up data. For the smallest problem size, the performance degrades pretty quickly. So if you look at the speed-up, it actually decreases when number of cores goes beyond 1,000. It is because the overhead of MPI operations quickly dominate the overall numerical cost when the local domain size is too small. But if we increase the problem size, then performance stays pretty good. It's very close to ideal value up to core number of 4,096. In terms of efficiency, it stays above 70%. Now we've covered the main features of SNAC. I'll give you a demo simulation to actually see how it works. So the example I brought here is oblique rifting problem. This is a model setting from Sandbox analog model. So on box filled with sand is considered as a block of crust, brittle crust, and it's pulled on two sides. And inside the domain, we have a rift zone, which has some obliquity to the extensional direction. And inside the rift zone, it can be strained, but outside of it, all the box remains as moves as rigid blocks. And on the right, you can see the surface deformation pattern after some extension. To set up a comparable model with SNAC, but in a real continental rifting scale, I set up a domain with the size of 100 by 100 kilometers and 20 kilometers thick. So this is a brittle crust block, and I'm pulling this block on two sides. As in Sandbox model, I define a rifting zone which can be strained, but outside of it, there's no gradient of velocity. And at the same time, I put some weak seed only within the rifting zone so that faults can develop out of those seeds. So this image shows the pattern of accumulated plastic strain after some extension. This is the animation that shows how actually the fault system develops. So you can see initial plastic seeds, and as extension goes on, you can see that faults grow out of one of those seeds, and they try to connect to each other. So the bottom box shows accumulated plastic strain distribution, and the top warped surface is the exaggerated topography. So the maximum relief from the highest to lowest point is about one kilometer. So let me just run the animation again. So individual faults are making pretty much a perpendicular direction to the extension. However, the interconnected structures are followed the predefined rifting zone. In that sense, it's pretty much consistent with Sandbox model. So in summary, snag is a 3D parallel Lagrangian code for elasto-viscule plastic material. It has been benchmarked and is decently documented. It's an open source being distributed through the website of computational infrastructure for geodynamics. And it's been linked to CSDMS website too. Now let's move on to a more important issue of how tectonic and surface processes interact with each other. So erosion in a convergent origin is now an iconic problem demonstrating those two processes interact in a non-trivial way. So this picture is from a paper by Willet, Giger 99, and it shows that pro or withdrawal wedge erosion can lead to completely different patterns of exclamation. So again, assuming that all of you are familiar with this type of interaction, I chose a couple of projects in the opposite setting. Opposite in a sense that we'll be looking at extension of tectonics rather than compression of one, and then sedimentation rather than erosion will be the process that adds complexity. So first example is the rifting style with and without sedimentation. It has been suggested that hot and thick content across when it's extended, then it will go through a wide rifting. Wide rifting meaning distribution of extension is pretty much uniform over a wide region. However, if the extended region is covered by sediments, then the style of rifting changes to narrow or focused rifting. Bialas and Bach showed this using a numerical model. So they set up a thick continental crust, which is about 50 kilometers thick and about three times wider. And if you look at the temperature field, it is hotter than the surroundings because it has a more abundant radiogenic heat. And as expected, the numerical model shows that without sedimentation, the rifting system evolves with a wide rifting style. So after about five million years of extension, which corresponded to 100 kilometers of extension, you can see from this phase field and temperature field that extended region has pretty much uniform thickness. And at the same time, if you look at the strain weight field, the two ends, two boundaries of the extending region is always the site of maximum deformation. However, you get completely different response when you start adding sediments. So on top of the extending region, this green stuff represents some sediment layer. So after about, again, 100 kilometers of extension, if you look at the strain weight field, now there is only one side of active deformation. And if you look at the temperature field, it has clearly asymmetric thickness variation. In other words, the whole system is now going through focused rifting. And the reason is that the thick blanket of sediments increased boundary layer thickness and the system will try to take advantage of non-uniform strength. In other words, the weakest part will only the weakest part survive and takes up all the extension. And another example I'm going to show is about rider blocks. Rider blocks is the pieces cut off from a hanging wall of a large offset normal fold. So it's a little tricky concept just described by words. So I'm going to show an animation. This animation captures the process of rider block formation. So this is 10 kilometer thick crust and 100 kilometer wide and is being pulled on two sides. And being supported by invested fluid on the bottom. So let me run it again. Normal fold develops and it rotates quickly. And as it rotates, it locks and then new fold forms. And at the same time, rider blocks form just one more time. So these blocks are called rider blocks. They are passively riding the foot wall and being transported. So here is a generic mechanism for rider block formation. So first you have new fold that follows so-called optimal angle. But as offset on the fold increases, the crust layers isostatically bend and therefore fold rotates and eventually it locks up. And if it cannot take any more offset, then the crust will develop a new fold at a higher angle than currently locked fold. Then a piece of crust that is between these two folds becomes a rider block. According to this conceptual model, we can infer two conditions for rider block formation. First we can say that a fold is too weak, then the fold is never going to lock up even if it becomes a low angle. And at the same time, if fold is strong enough if there is nothing on the hanging wall side, in other words, if there is no sedimentary infill in this rifted valley, then rider blocks wouldn't form. So if we turn these arguments around, then the presence of rider blocks implies that fold is strong enough and there should be enough amount of sedimentary infill. And these two conditions for rider block formation has very interesting implication. So I'm now showing you oceanic core complex like its continental counterpart in an oceanic core complex, mid-crustal material is exposed on the surface through a large ups and normal fold. So this specific core complex is at about 30 degrees north along mid-Atlantic ridge. And this is an Atlantis transform and it's right north of it. So it's called Atlantis massive. So if you look at the more detailed bathymetry map, this whole block is called Atlantis massive. And if you take cross section of topography and also structure interpretation from south to north, then it looks like this. So at the southern tip, it has a very smooth and corrugated surface. But to the north, it has more rugged topography. The interpretation of this difference, a long-rich difference in topography by restaurant Ranero is that as you go north, you get closer to segment, rich segment center where magnetism is more active. That means the whole rift valley is filled with volcanic infill but at the tip of rich segment, magnetism is not as active, so it doesn't have that much of volcanic infill. That means even if fault can be strong enough to have right of blocks, at the southern tip where there is no sedimentary infill, you cannot have right of blocks. But if you get closer to segment center, then there is enough of sedimentary infill and fault is strong enough so that you can have multiple right of blocks. So we've just seen that sedimentation or space filling surface process is also important in tectonic processes. So in 2D settings, sediments are usually assumed that it's transported along the third dimension in and out of plane and as much as we need. So a lot of hand waving is going on there. And as we've seen in the oceanic coral complex problem, some problems are inherently three dimensional. So that leads us to conclude that physically based 3D coupling between tectonic and surface process modeling is desired. There is one ongoing effort led by Federal Opton in JNS Science, New Zealand and Greg Tucker here at University of Colorado. They are trying to couple a tectonic modeling code, the FLAC3D, that is a commercial geotechnical modeling code and child, the surface processing code. So as a pilot model, they set up a 60 by 60 kilometers and 10 kilometers thick domain, which is pushed from the left at a rate of 1 centimeter per year and the whole left half of the domain is actually kinematically constrained and the right half has a frist lip boundary condition. So according to these kinematic boundary conditions, mountain belt will be built in the middle of domain and that topography developed in the code FLAC3D is given to child and child will adjust that topography by erosion and that adjusted topography is given back to FLAC3D and keep going that way. So here is a preliminary result from this coupling project. On the left, the material is assumed highly erodable and on the right, material is less erodable. The result you are supposed to be seeing is that highly erodable material has slightly larger wavelength of topography along the mountain belt and less erodable material has higher frequency topography and at the same time the maximum topography for less erodable material is obviously a lot higher. So this simplistic and preliminary model already shows that this is the power of code coupling. Then the question is can SNAC be used in place of FLAC3D in this type of project? So FLAC3D's 3D SNAC is 3D2. SNAC can handle inelastic deformation which is necessary for build a mountain belt and SNAC can handle non-uniform boundary conditions too. The remaining question is can SNAC be coupled with child? The answer is actually yes but the real question should be asked is how? And SNAC is based on software framework so I believe it's already in a pretty good shape to be incorporated into CMT system. So I'd like to come up with a little more concrete plan for this coupling during this meeting and use of open source like SNAC should be beneficial not only financially but also scientifically because researchers can modify the code for their own purposes. Thank you for your attention. The seeding that you did for the first example of SNAC, how much does that determine the outcome of the ripped geometry, number of rifts, spacing and so on? It's very hard to quantify. So I'm actually trying to calculate collective property like the displacement versus, well the length versus the fault drop ratio or the size distribution. So it can be defined only that in a collective way. So if I change the distribution of initial seeds it will definitely change the pattern but it's not clear whether that collective property will be also changing. Any other questions? Great. Can you tell us a little bit about whether SNAC can go to very high strains? Do you use re-gridding or something like that? Yes, I do, yes. Well to add something there is issue with re-gridding. As you know re-gridding involves numerical diffusion which can be sometimes very serious. So it would be great if you can come up with the better algorithm than the currently implemented SNAC. Again a simple question but when you actually put the sediment, the example that you showed where you filled your basin with sediment, where did that sediment come from? Do you have like rudimentary surface processes that you actually removed material from the block and deposited or you just basically created mass? Short answer is we simply create material as much as we want. But that model actually included very simple like the diffusion style topographic modification but that amount is not totally constrained, not totally constraining the amount of actual sediment we put in the basin. Thank you very much.