Category Archives: geoscience

Building a simple delta numerical model: Part III

In this part of “Building a simple delta numerical model”, we’ll simply develop a module for calculating sediment transport at all locations within our model domain. This is naturally a critical piece of the model to get right, because the magnitude of sediment transport that occurs within the delta system will have a direct correlation with the magnitude of geomorphic change to the delta system. Just as before, we’ll be working our module for the Mississippi River, and so we will use the Engelund and Hansen, 1967 formulation for bed-material load sediment transport.

It is worth noting here that there are countless sediment transport equations published within the literature, with applicability for different scenarios, covering the entire range of alluvial rivers, the shelf environment, to turbidity currents. A nice review of some relations is given by Garcia, 2008, who produces the plot below with six selected relations, where τ* is shear stress (Shields number) and q* is a dimensionless sediment transport.

garcia_formulas

The only points I hope to make with this figure are that stress and transport are nonlinearly related, and that many of the relations are quite similar. In fact, the majority of sediment transport equations are based on an excess-shear stress formulation, where, as shown above, the sediment transport prediction scales with the shear stress. The Engelund and Hansen, 1967 equation is given in its dimensional form by (Equation 1)

eh_formula

where qs is the sediment transport per unit width (m2/s), β is an adjustment coefficient discussed below, R = 1.65 is submerged specific gravity of the bed sediment, = 9.81 m/s2 is gravitational acceleration, is the grain size in question, Cf is the coefficient of friction, ρ = 1000 kg/m3 is the fluid density, and τ = ρ Cf U2 is the fluid shear stress, where is depth-averaged fluid velocity.

The solution to this equation is obtained by simple algebra, and so I will display the code calculations below without any explanation.


U = Qw ./ (H .* B0); % velocity

function [qs] = get_transport(U, Cf, d50, Beta)
% return sediment transport at capacity per unit width
% 1000 = rho = fluid density
% 1.65 = R = submerged specific gravity
% 9.81 = g = gravitation acceleration
u_a = sqrt(Cf .* (U.^2));
tau = 1000 .* (u_a.^2);
qs = Beta .* sqrt(1.65 * 9.81 * d50) * d50 * (0.05 / Cf) .* ...

        (tau ./ (1000 * 1.65 * 9.81 * d50)) .^ 2.5;
end

Let’s first simply examine the behavior of Equation 1, in its dimensionalized form, for the Mississippi River. We need to define the grain size and adjustment coefficient to use in our model; we’ll follow from the literature and use 300 μm for the grain diameter, assuming this makes up the entirety of the bed, and β = 0.64 an empirical adjustment (Nittrouer et al., 2012; Nittrouer and Viparelli, 2014).

Evaluating the equation over a range of velocities from 0 to 3 m/s, we find that the nonlinear behavior shown in the above figure is also true for the Engelund and Hansen formulation.

 

Now, let us implement this into our delta model. Following from earlier, the flow depth and velocity are related by the cross sectional area and discharge, where for a rectangular cross section, the velocity = Qw / H B. We therefore can take our calculation of flow depth from the previous lesson, and solve for velocity at each location. We then simply plug this into the relation for τ given above and solve Equation 1.

 

The bottom plot reveals the behavior of the flow velocity calculated through the backwater region (shown in red), and the resulting sediment transport calculated (shown in green). The nonlinearity between stress (i.e., velocity) and transport is again revealed by the change in slope of the lines through the backwater region. Following from left to right along the green curve, it becomes obvious that there is a decrease in transport downstream. If there is more sediment being transported upstream than downstream, then there must be sediment that is deposited to the bed somewhere along the channel bed between ‘upstream’ and ‘downstream’.

This thought experiment represents one half of the famous Exner equation, which we will explore in the next part, and will provide the last module to complete our delta model. My complete code for the sediment transport calculation in our backwater delta model can be found here.

 

References

1. Engelund, F. & Hansen, E. A monograph on sediment transport in alluvial streams. (Technisk Vorlag, Copenhagen, Denmark, 1967).
2. García, M. H. Sediment Transport and Morphodynamics in Sedimentation engineering processes, measurements, modeling, and practice 21–163 (American Society of Civil Engineers, 2008).
3. Nittrouer et al., Spatial and temporal trends for water-flow velocity and bed-material sediment transport in the lower Mississippi River. 2012. Geological Society of America Bulletin, 5 (Table 1).
4. Nittrouer, J. A. & Viparelli, E. Sand as a stable and sustainable resource for nourishing the Mississippi River delta. Nature Geoscience, 7, 350–354 (2014).

This material is based upon work supported by the National Science Foundation (NSF) Graduate Research Fellowship under Grant No.145068 and NSF EAR-1427177. Any opinion, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

Building a simple delta numerical model: Part II

Taking the simpler form of the backwater equation, that is, assuming that width does not vary over the channel reach, we find an expression for changing flow depth along the channel streamwise coordinate x: (Equation 1)backwater_fixed_eqnwhere S is the channel bed slope, Cf  is the dimensionless coefficient of friction, and Fr is the Froude number: (Equation 2)

A couple quick notes about the overall delta model we are developing herein:

  1. The code for this entire model will be developed in Matlab syntax, however, the math is all just math, and the code could easily be translated to another programming language.
  2. The model will be set up for the conditions of the Mississippi River, to what I feel are the best constraints published in the literature. Surely not everyone will agree with some of my choices.

We set up our model domain by defining some parameters:

L = 1200e3; % length of domain (m)
nx = 400; % number of nodes (+1)
dx = L/nx; % length of each cell (m)
x = 0:dx:L; % define x-coordinates
start = 63; % pin-point to start eta from (m)
S0 = 7e-5; % bed slope
eta = linspace(start, start - S0*(nx*dx), nx+1); % channel bed values
H0 = 0; % fixed base level (m)
Cf = 0.0047; % friction coeff

bw_schem_label

Each of these terms are necessary for defining the boundary conditions the backwater calculation requires. We will additionally need to define the water discharge (Qw), and it’s width averaged derivative qw = Qw / B, where B is the channel width. We’ll use a fixed width of 1100 m (Nittrouer et al., 2012) and a discharge of 10,000 m3/s; since width is fixed in our model, qw = 31.8182 m2/s.

Finally, to solve the above backwater equation, we’ll need to find the local slope at all nodes — our bed holds a constant slope so we could just skip this step and send S0 directly into the calculation, but we’ll be wise beyond our present experience, and prepare our calculation now for the future when our bed slope is constantly changing in time and space. A simple slope calculation is completed by a central difference calculation, applying forward difference and backward difference calculations at the domain upstream and downstream boundaries, respectively.

function [S] = get_slope(eta, nx, dx)
% return slope of input bed (eta)
S = zeros(1,nx+1);
S(1) = (eta(1) - eta(2)) / dx;
S(2:nx) = (eta(1:nx-1) - eta(3:nx+1)) ./ (2*dx);
S(nx+1) = (eta(nx) - eta(nx+1)) / dx;
end

As described in the previous section, we will use a predictor-corrector scheme to find a solution to Equation 1. This is not meant to be a comprehensive guide on solving ODEs, and so I will only summarize the method here.

  1. We first can evaluate the flow depth H at the downstream boundary, where elevation is fixed (H = H0 – η = 21 m).
  2. Then we predict the value at the next upstream node by solving Equation 1 for dH/dx at our known downstream-most node, using the inputs described above and H and S at the known node. This yields dH/dxp = 6.5784e-05.
  3. Multiplying dH/dx by dx = 3000 m and subtracting this change from the known H at the downstream-most node, we find our predicted value Hp for the next upstream node.
  4. We now repeat the evaluation of Equation 1 with the same inputs except using our new Hp value in Equation 2. This yields our corrected change in flow depth over x:
    d
    H/dxc = 6.5663e-05.
  5. We then combine the results of these two predictions in flow depth change to give our best solution to the ODE Equation 1 at the downstream node. We combine the values using the trapezoidal rule for Hi = H(i+1) – ((0.5) * (dH/dxp + dH/dxc) * dx), giving a predicted change in flow depth of 0.1972 m, for a depth at Hi of 20.8028 m (it was 21 m at the downstream node).
  6. repeat steps 2-5 solving for each node moving up the model domain until the final node is determined.

This method is fairly stable for most discharges, but may fail when discharge is very low. In thsi case, the backwater region becomes condensed, and the rate of change in H becomes very large over a very small distance. My code for the predictor corrector implementation can be found as the get_backwater_fixed() function within the complete backwater_fixed model. Below is a gif showing the results of a range of input discharges to our model.

 

and a still multi-surface example to demonstrate the variability.

 

The low and moderate discharge water surface curves in the above figure (solid and dashed lines) demonstrate the classical convex-up (aka M1, e.g., Chow, 1959) condition that characterizes lowland delta systems during most of the time. The third, high discharge, water surface profile exemplifies a convex-down water surface profile. This condition, called the M2 condition has been recently demonstrated to have importance to the cyclic processes of avulsion and delta-lobe growth (e.g., Chatanantavet et al, 2012; Ganti et al., 2016).

 

Building a simple delta numerical model: Part III

 

References:
1. Nittrouer et al., Spatial and temporal trends for water-flow velocity and bed-material sediment transport in the lower Mississippi River. 2012. Geological Society of America Bulletin (Table 1).
2. Chow, Ven Te. Open Channel Hydaulics. N.p., 1959. Print. McGraw-Hill Civil Engineering Series.
3. Chatanantavet, Phairot, Michael P. Lamb, and Jeffrey A. Nittrouer. Backwater Controls of Avulsion Location on Deltas. 2012. Geophysical Research Letters 39.1
4. Ganti, V., A. J. Chadwick, H. J. Hassenruck-Gudipati, B. M. Fuller, and M. P. Lamb. Experimental River Delta Size Set by Multiple Floods and Backwater Hydrodynamics. 2016. Science Advances 2 (5).

This material is based upon work supported by the National Science Foundation (NSF) Graduate Research Fellowship under Grant No.145068 and NSF EAR-1427177. Any opinion, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

Building a simple delta numerical model: Part I

This is the first in a series I’m going to write over an undetermined-number-of-parts series where we will develop all the pieces to a simple delta numerical model. I’ll outline all the pieces here while I work up an interactive version with a GUI in Matlab, and we write pieces of the code together within the posts.

——————————————————————————–

 

The water depth (H) over a channel bed with a fixed-elevation downstream boundary can be calculated by a mass-and-momentum conservative derivation following from the full Navier-Stokes equations (e.g., Chow 1959). The so-called backwater equation, in one form, is presented as (Equation 1)
backwater_fixed_eqnwhere H is the flow depth, x is the streamwise coordinate of the flow, S is the channel bed slope, Cf  is the dimensionless coefficient of friction (corresponding to form drag in the channel bed and walls (i.e., roughness and bedforms)), and Fr is the Froude number, which scales with the ratio of inertia of the flow to the gravitational field driving the flow and can be given as (Equation 2)where qw = Qw is the width averaged water discharge, and g is the gravitational constant. Equations 1 and 2 are derived for a width-averaged, constant flow width scenario.

Equation 1 is useful to the research we do as scientists of lowland depositional systems because it predicts changes in flow depth through a delta system. The equation is solved with a boundary condition of a known flow depth (H) at the fixed-elevation downstream “receiving basin” boundary. For variations in river conditions (i.e., river water discharge), the surface elevation of receiving basins of large lowland rivers are effectively fixed, therefore the equation can be used to estimate a water surface over the lower reaches of the rivers we study for any input discharge conditions!

bw_schem_simp

With a calculation of H, we can then use a conservation of water equation Qw = U A (where U is flow velocity and A = B H is channel cross sectional area for an approximately rectangular channel), to evaluate the flow velocity at all points along the channel. Calculation of U in turn allows for calculation of sediment transport, and ultimately sediment deposition — but more on these parts in future parts of this series!

Equation 1 is valid for a solution where width (B) is fixed over the downstream length of x. Where width is roughly constant, say within the lower reaches of a channel we can use this simple form of the backwater equation. However, where width is not constant over all x (i.e., dB/dx != 0), say in the river plume beyond the channel mouth, we must retain an addition term from the derivation from the Navier-Stokes equations to account for this change in width. Since we’re using computers to solve the equations, and the dB/dx term is basically only one extra gradient calculation, this will an easy addition to our model in the future. I’ll write a separate post and link it here with a solution for that version of the backwater equation as well.

So how do we solve Equation 1? Let’s see what we know: the elevation of the receiving basin is fixed at z = 0. Therefore, if we know the channel bed profile, we know at the downstream boundary for ALL possible values of Qw. We can evaluate Equation 1 from the known downstream boundary, and knowing all variables on the right hand side at this boundary, calculate the change in flow depth moving upstream over the channel bed (dH/dx)!

The solution to Equation 1 (or in its slightly more complex form for a varying width domain) can be numerically estimated pretty easily by a predictor-corrector scheme. We’ll work up a calculation for a fixed-width form of the backwater equation (Equation 1) in the next post, and then proceed to use this as one module of our simple delta model.

Building a simple delta numerical model: Part II

 

This material is based upon work supported by the National Science Foundation (NSF) Graduate Research Fellowship under Grant No.145068 and NSF EAR-1427177. Any opinion, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

Dietrich Settling Velocity Matlab code

Just about two years ago, I published a code for calculating the settling velocity of a particle in water by the Dietrich, 1982 method. I later realized an error in the code that made it produce erroneous estimates. I replaced it with a corrected version of the code written in R, that you can still find on my Github page, but won’t link here. I’ve taken down the old Matlab code and unpublished the original post.

Since the first publication two years ago, I’ve dramatically improved my technical coding skills, and, just as importantly, learned a ton about how to share code with others. A new version of the code can be found linked a the end of this post, but first some philosophy about coding.

As scientists, we want to be able to reliably reproduce results and we have an attraction to openness and access to knowledge, as such, we inevitably share code with one another. I’ve shared a good bit of my code thus far in my career, and I’ve made a few notes about the matter below:

  1. It is better to organize your code in a way that may not be the most succinct (when speed is not an issue) while improving the readability of the code. In essence, making your code readable to those unfamiliar with it, so they can follow along with your logic.
  2. Comment everything. It is so easy to get into the groove of writing something and blow through it until it’s finished and then never return to clean up and document the code. When others (or even you) go back to look at it, simple comments explaining what a line does go miles further than a detailed explanation in the documentation.
  3. Write your code to withstand more that you are currently using it for. What I mean by this is to not write code that works only when you are manipulating this particular dataset or working this particular problem, but code that will work over a broader range of data or can be applied to another problem. Along the same line, use explicit values as sparingly as possible, e.g., if X = [1:12] is a 12 element vector and you want the middle value, do not use X(6), but instead use X(length(X)/2); or slightly better yet, X(ceil(length(X)/2)) to at least not throw an error and give you the first of the two middle rows if X is an odd number of elements in length.

 

Anyway, the code is uploaded to my Github and you can find it there. The documentation output on a help get_DSV call in Matlab is produced below.

get_dsv_documentation

 

References:
Dietrich, W. E. Settling velocity of natural particles. Water Resources Research 18, 1615–1626 (1982).

This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No.1450681. Any opinion, findings, and conclusions or recommendations expressed in this material are those of the authors(s) and do not necessarily reflect the views of the National Science Foundation.

Groundwater flow demonstration slides

I recently came across a set of old slides from the National Water Well Association, created in 1977. The slides can be found here as an imgur album. These images show principles of groundwater flow that were not novel at the time, but the images do a great job of showing (through time series primarily) how different mediums can drastically change groundwater flow. Further, the slides show the effects of groundwater extraction through natural (gaining streams) and anthropogenic (wells). Below is a gif of the first ~20 slides I put together, which shows the experimenters setup over the duration of the experiment. There is a constant “head”, or a presence of water, at the right side that is flowing towards the left side of the experimental setup, with some dye being injected at a few depths to show the flow lines.

groundwaterexperiment

I’ve included my favorite image from the set below, which shows (beautifully) a cone of depression. When initially drilling a well, it (obviously) is critical to put the well into the water table, such that the pumping end of the well is entirely submerged. Once the well is “switched on”, there is a drawdown of the water table above the well pump head.

groundwater14

The cone forms around the well extraction point with its deepest point at the well head, as the well pulls water in from ALL directions (below included!). This disturbs the flow lines, which is shown really well by this experiment. Note that this experiment just shows a 2D slice of the cone, when in reality the “aquifer” is 3D and would extend into and out of your computer screen, making the wedge you see here a cone in three dimensions.

Yellow River shoreline traces

For some ongoing research I needed a set of shoreline traces of the Yellow River delta that I could measure some properties of to compare modeling results to. My model predicts delta lobe growth rates, for example the growth of the Qingshuiguo lobe on the eastern side of the delta from 1976 to 1996.

early growth stage of Qingshuiguo lobe of the Yellow River delta, figure from van Gelder et al., 1994.

early growth stage of Qingshuiguo lobe of the Yellow River delta, figure from van Gelder et al., 1994.

It’s a little early for me to release the data I’ve collected from the shoreline traces (n=296), or the code I used to generate the traces — but below is a figure mapping all the shoreline traces I have from 1855 to present. My model can also predict rates of delta system growth (i.e., growth over many lobe cycles) and I have generated data from these traces that I can compare to as well. I like the below figure, because I think it nicely demonstrates the rapid reworking of delta lobes (yes even on “river dominated” deltas!); evidenced most clearly by the retreat of the Diakou lobe to the north side of the delta.

xom_shoreline_map

 

This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No.1450681. Any opinion, findings, and conclusions or recommendations expressed in this material are those of the authors(s) and do not necessarily reflect the views of the National Science Foundation.

Yellow River (Huanghe) 2016 survey plans

I’ll be departing soon for my second field campaign on the Yellow River, China. We traveled there around the same time last year to generate a dataset we could begin to explore and develop research ideas; now, a year later, we will return with a more focused plan to gather the data we need to address our research questions. We will bring with us a suite of sophisticated data collection equipment. Below is a terrible sketch I made to demonstrate most of the data we will collect at each of our stations located on the river.

YR2016_deployment

Pretend that the grey thing floating on the water might look like a boat, viewed from the back. We aim to characterize the mechanisms for sediment transport within the Yellow River at various stages of a flood discharge curve. We use a point-integrating sediment sampler to collect numerous suspended sediment and water samples from various depths in the water column (yellow stars). In order to compare the compositions of suspended material to its source (i.e., the bed), we use a Shipek grab sampling device to collect sediment from the bed for our analysis. Deployed off the other side of the boat, we will can characterize the velocity structure of the water column through the use of an aDcp (acoustic velocity profiler) and a mechanical propeller driven velocimeter for near-bed measurements where the aDcp may lack resolution.

This setup represents only a portion of our survey plans, of course we’ll be there for six weeks that will generate a diverse and (hopefully!) comprehensive dataset for our future work.

Hiking Garmisch-Partenkirchen

I’m posting this in the off-chance it may help some future traveler who wants to find the same kind of day hike I did, so that he/she can enjoy it as much as I did. I recently spent a few days in München and wanted to take a day trip to the Alps for a hike. I was having an awful time finding any actual hiking routes, so I just kind of figured I would wing it once I got there. So I boarded a train the south for a city called Garmisch-Partenkirchen. On the train I asked a fellow rider (German) if they knew of any hiking trail out of G-P. While not immensely helpful, he did suggest to me my starting point.

2016-05-14 16.25.12

The train station at Garmisch-Partenkirchen

The idea was to set out from the train station in G-P for a gorge called Partnachklamm (the Partnach river is the river that has cut the gorge); from there I would just continue climbing and wandering up trails and the mountains until I got tired and then head back down. The Partnachklamm begins about 1 km up a small road that begins at the Olympic Ski Jump facility. The hike I followed was ~19 km, or approximately 10 km up and 9km down. HERE is a .kmz file so that you can see the approximate path I followed, you may have to click “view raw” in order to initiate the download (please note that I just mapped this now in Google Earth, and it does not necessarily map the exact trails I followed!!).

So for some direction, the hike I enjoyed could be mapped as follows (be sure to enjoy the scenery along the way too!).

  1. Head out to the north and east of the train station, along St. Martin-Straße toward the east. When you near the river, head south on Partnachhauenstraße which runs along the river. Here there is a neat “geology of the Alps” exhibit that actually has some spectacular specimens. The path crosses the river at Schornstraße. Continue to follow the path to the ski jump stadium.
  2. To the west of the stadium (right hand side as seen from your approach) is a paved road. Follow it (and signs) for the Partnachklamm. There is a 4 euro entrance fee to the Partnachklamm, which may seem unfair, but it was definitely worth while to me. Enjoy the hike through the Partnachklamm.
  3. Once out of the gorge, you will follow the path up a series of switchbacks, following signs for “Das Graseck”. This will take ~30 minutes.
  4. Continue beyond Das Graseck and take the paved road up the hill on your right. Eventually the pavement ends in favor of a dirt path at the “peak” of a hill. You then want to follow the small dirt path up, in the steepest direction up a long series of switchbacks through a dense pine forest. This will be the path to the “Berggasthof Eckbauer”. This will take ~45 to 60 minutes.
  5. Finally you will reach the Berggasthof Eckbauer. Have a pint and enjoy the spectacular view.2016-05-14 14.03.29
    2016-05-14 14.09.20
  6. My final piece of advice is when you return down, you will want to follow the same path you followed up, but once you pass Das Greseck, you will follow a steep path down and to the right. This path will lead you down to some waterfalls and you will see a sign warning that the trail is for experienced hikers only. You should heed this warning, but also don’t be frightened by it. This path will return you to the “tollbooth” of the Partnachklamm but allow you to walk above the gorge and get some different perspective.

Lastly, here is my new computer background. Go and hike the spectacular German Alps.

2016-05-14 15.01.11

Gastroliths in the Morrison Fm

I’ve recently returned from a trip to New Mexico for the field methods class of which I am a TA. It was an awesome trip and great experience for me in teaching, but that’s another story. The field area we work in is within the Jurrassic Morrison depositional basin (active roughly during the Kimmeridigan ~157-152 Ma). Within the Morrison Formation is the Brushy Basin member (abbreviated Jmb), renowned for the abundance of dinosaur fossils found within the rock unit. Jmb is found within our field area, so I told the students to keep an eye out for any good finds when walking with the unit.

Although Jmb has an abundance of fossils, we didn’t find any. BUT, the Morrison is also famous for its bounty of another paleontological tool, the gastrolith. Gastroliths are interpreted as stones that dinosaurs would ingest in order to aid with breaking down food and aiding in digestion. A gastrolith may be more generally defined as “a hard object of no caloric value (e.g., a stone, natural or pathological concretion) which is, or was, retained in the digestive tract of an animal” [1].

There are gastroliths all over within the Jmb, some small, and some larger. Below are three photos of the “best” gastrolith I found on our trip.

gastrolith3 gastrolith2 gastrolith1

Gastroliths are often recognized by their very smoothed and polished appearance (some other examples here). I suppose that to be certain my rock is a gastrolith, and not simply a rock polished by water or wind, it should be found in association with the remains of the animal it was within. Regardless, I’m really glad to be able to add a rock with such an interesting back-story to my collection.

 

 

Adding real data into GMT map, making a hillshade, and positioning

Alright, let’s get down to business plotting some real data on a map. Grab the zip file here, and uncompress it in the folder you want to build a map in. The zip contains the data to map, as well as  on of my favorite color palettes from cptcity by M. Burak Yilkilmaz (file named mby.cpt).

Here’s the end product of this tutorial, we’re going to make two maps, the first will be a simple data map, and the second a hillshaded version of the same map. We’ll also learn how to position them side by side like below.

g12

We’ll start by plotting a simple map, with just the data represented as different colors (left side of the map). The data we’re plotting are from a digital elevation model (DEM), and of Eastern North America. Making this image is simple, we simply need to call
gmt grdimage real_data.nc -R278/286/36/42 -JM4 -B2WESN -Xc -Yc -Cmby.cpt -V > flat.ps
where you already know most of the switches above. grdimage is a gmt tool that you can use to plot any data stored in the netCDF format. To color the image based on the values in the .nc file, we need to use a color palette. The color palette is called with the -C switch, followed by the name of the color palette. More on color palettes in a future tutorial though.

Now, let’s talk about hillshades. A hillshade (also known as shaded relief) is an element of a map that adds depth to topography. It is essentially simulating a light source and the shadows cast on the landscape from the hills onto the valleys. Making one in GMT is straightforward but requires a few steps. The main step of the process is done first with grdgradient. We use the following command
gmt grdgradient real_data.nc -Ghillshade-grad.nc -A345 -Ne0.6 -V
In this command, -G has a different meaning than in this tutorial; another common use of the -G switch is to give the name of the output file for the command. -A is the azimuth from which the light source should be directed, and -Ne is a normalization parameter for the output from the calculation, just use 0.6.

Then, we will ensure that the shaded relief map is reasonably balanced by normalizing it again along a Gaussian distribution
gmt grdhisteq hillshade-grad.nc -Ghillshade-hist.nc -N -V

Finally, in order to ensure the values are between -1 and 1, a final step required for creating the intensity file needed for a hillshade to the image. We do this by dividing by some value just larger than the range of values in the histogram equalized grid file. Let’s find this range by looking at the info for this file.

$ gmt grdinfo hillshade-hist.nc
hillshade-hist.nc: Title: Produced by grdhisteq
hillshade-hist.nc: Command: grdhisteq hillshade-grad.nc -Ghillshade-hist.nc -N -V
hillshade-hist.nc: Remark: Normalized directional derivative(s)
hillshade-hist.nc: Gridline node registration used [Geographic grid]
hillshade-hist.nc: Grid file format: nf = GMT netCDF format (32-bit float), COARDS, CF-1.5
hillshade-hist.nc: x_min: 278.004166667 x_max: 286.004166667 x_inc: 0.00833333333333 name: longitude [degrees_east] nx: 961
hillshade-hist.nc: y_min: 35.9958333333 y_max: 41.9958333333 y_inc: 0.00833333333333 name: latitude [degrees_north] ny: 721
hillshade-hist.nc: z_min: -4.67873573303 z_max: 4.67873573303 name: Elevation relative to sea level [m]
hillshade-hist.nc: scale_factor: 1 add_offset: 0
hillshade-hist.nc: format: netCDF-4 chunk_size: 138,145 shuffle: on deflation_level: 3

So we will use 5 as our divisor, in the command
gmt grdmath hillshade-hist.nc 5 DIV = hillshade-int.nc
which takes our -hist file and divides every grid cell by 5, putting it within the range of -1 to 1.

Finally, we use
gmt grdimage real_data.nc -R278/286/36/42 -JM4 -B2WESN -Xc -Yc -Ihillshade-int.nc -Cmby.cpt -V > shaded.ps
the grdimage tool again, in almost exactly the same way as before, except adding the -I switch with our normalized intensity file as the input.

Now, to put the two images side by side, we need to set different -X and -Y targets for each image. These switches set the offset of the lower-left corner of the PostScript layer being produced, from the last reference point used on the PostScript page. The a argument can be added to the switches to reset the reference point to the lower-left corner of the page. Try using -Xa1 -Yc for the flat image and then overlay the hillshade image at -Xa6.5 -Yc. Remember, you’ll need to add in -K, -O, and >> in order to produce a multi-layer map.

See if you can make it work by yourself, but if you need a hint, you can get the script I used from here.