How Often Does Nepal Shake?

Start with Since: 1973 and Chart: Ordinary.

• The ordinary chart. The count of magnitude 4 or more is 1,511; of 6 or more, 17. On an ordinary scale the small earthquakes tower over the rest, and everything from 6 up lies flat along the bottom, as if nothing happened there.
• Log scale. Now each line up the chart is ten times the one below: 1, 10, 100, 1,000. From magnitude 4.5 to about 6.5 the dots fall close to a straight line, and the dashed line beside them goes down by ten times for each step of one magnitude.
• Move Magnitude from 5.0 to 6.0 to 7.0, and read How often: about 4 a year, then one every 3 years, then one every 18. The map’s red rings show where they were: nearly all along the Himalaya, in a band through Nepal from west to east.
• Since: 1900. The great earthquake of 1934 comes back, and the count of 6 or more doubles, from 17 to 35. But the count of 4 or more hardly grows, from 1,511 to 1,578: 73 more years added only 67 earthquakes. Look at Since: 2000 too.
• Magnitude 4.0. On the log scale, the dots from 4.0 to 4.4 bend away below the dashed line.

What a magnitude is

A magnitude is worked out from the shaking that seismometers record, allowing for how far away they are. It goes up in steps that multiply: one step up is about 10 times the ground’s movement and about 32 times the energy released (as the USGS explains it). So a magnitude 7 releases about 32 × 32, roughly 1,000 times, the energy of a magnitude 5.

Magnitudes come out a little different from different measurements and agencies. The earthquake of 25 April 2015 was 7.8 by the USGS, and 7.6 as Nepal’s National Seismological Centre (NSC) recorded it; the great earthquake of 1934 is 8.0 in the USGS’s catalog and 8.4 in the government’s report after 2015. The file in this topic has the USGS’s, so its numbers agree with each other.

Counting “this big or bigger”

You could count how many earthquakes were exactly of each magnitude: 4.5, 4.6, 4.7… But with only a few big earthquakes, those counts jump about: one magnitude might have 3, the next 0, the next 2. Seismologists count instead how many were of each magnitude or more. Each count then includes all the bigger ones, so it can only go down as the magnitude goes up, and it goes down smoothly.

A straight line on a log scale

When a count gets ten times smaller at every step, an ordinary chart squashes everything but the first step against the bottom. On a log scale, each step of ten times is the same height, so a count that is divided by ten at every step goes down in a straight line. A log scale says it with the logarithm: log₁₀ of 1,000 is 3, of 100 is 2, of 10 is 1, of 1 is 0. Dividing by ten takes one off the logarithm.

So the law is a straight line between the logarithm of the count, N, and the magnitude, M:

$$\log_{10} N = a - bM$$

a says how many earthquakes there are in all (it depends on the place and on how many years are counted), and b, the b-value, is how steeply the line goes down. If b is 1, each step of one magnitude divides the count by 10¹, ten times. In most places with many earthquakes, b comes out close to 1.

The data is the file nepal_quakes.csv: one row per earthquake, its date, its magnitude and where it was. Open it to see it.

Step 1: the biggest

import pandas as pd quakes = pd.read_csv("nepal_quakes.csv")print(quakes.shape)print(quakes.head(3))print()print(quakes.sort_values("mag", ascending=False).head(6))

(1578, 4)         date  mag    lat    lon0  1911-10-14  6.6  30.99  80.341  1913-03-06  6.5  30.20  81.932  1916-08-28  7.0  29.67  80.68             date  mag    lat    lon9     1934-01-15  8.0  27.27  86.94867   2015-04-25  7.8  28.23  84.731021  2015-05-12  7.3  27.81  86.071435  2025-01-07  7.1  28.60  87.382     1916-08-28  7.0  29.67  80.6816    1936-05-27  6.9  28.38  83.32

1,578 earthquakes. The biggest is the great earthquake of 15 January 1934 in eastern Nepal; then Gorkha’s, 25 April 2015; its biggest aftershock, 12 May 2015, near Dolakha; and the earthquake of 7 January 2025 in Tibet, north of the Everest region. The numbers on the left are each row’s place in the file.

Step 2: which years can be counted?

The think box above said the old years miss the small earthquakes. Count them, decade by decade. date is text such as "2015-04-25", so .str[:4] takes its first four letters, the year.

import pandas as pd quakes = pd.read_csv("nepal_quakes.csv")quakes["year"] = quakes["date"].str[:4].astype(int)quakes["decade"] = quakes["year"] // 10 * 10 small = quakes[quakes["mag"] < 4.5]big = quakes[quakes["mag"] >= 5]table = pd.DataFrame({    "4.0 to 4.4": small.groupby("decade").size(),    "5 or more": big.groupby("decade").size(),}).fillna(0).astype(int)print(table)

        4.0 to 4.4  5 or moredecade                       1910             0          61920             0          21930             0         141940             0          21950             0         121960             0         231970             9         281980            27         211990           126         332000           163         372010           307         672020           173         34

The earthquakes of 5 or more come at much the same rate from the 1960s on (the 2010s have more: Gorkha’s aftershocks). The small ones of 4.0 to 4.4 appear only in the 1970s and keep growing for decades: the earth didn’t shake more, the world’s seismometers caught more. So from here on, the programs use the earthquakes since 1973, and fit the line only from magnitude 4.5, where the counts look complete.

Step 3: the counts, two ways

For each magnitude, count the earthquakes since 1973 of that magnitude or more. recent["mag"] >= m is a column of True and False, and .sum() counts the Trues. Then draw each count as a bar of #s, first by the count (one # for 25 earthquakes), then by its logarithm (math.log10, ten #s for each ten times).

import mathimport pandas as pd quakes = pd.read_csv("nepal_quakes.csv")quakes["year"] = quakes["date"].str[:4].astype(int)recent = quakes[quakes["year"] >= 1973] for m in [4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5]:    n = (recent["mag"] >= m).sum()    print(f"{m} or more: {n:5}  {'#' * round(n / 25)}")print()for m in [4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5]:    n = (recent["mag"] >= m).sum()    print(f"{m} or more: {n:5}  {'#' * round(10 * math.log10(n))}")

4.0 or more:  1511  ############################################################4.5 or more:   706  ############################5.0 or more:   212  ########5.5 or more:    46  ##6.0 or more:    17  #6.5 or more:     9  7.0 or more:     3  7.5 or more:     1   4.0 or more:  1511  ################################4.5 or more:   706  ############################5.0 or more:   212  #######################5.5 or more:    46  #################6.0 or more:    17  ############6.5 or more:     9  ##########7.0 or more:     3  #####7.5 or more:     1

In the first chart, the earthquakes from 6.5 up get no # at all. In the second, the bars get shorter by about the same amount at each step: the staircase of a straight line. The step from 4.0 to 4.5 is the shortest, the small earthquakes the catalog misses. (The bar for 7.5 is empty because log₁₀ of 1 is 0.)

Step 4: fit the line

A computer can find the straight line that comes closest to the points: for each point it measures how far the line misses it, up or down, squares the miss, and chooses the line that makes those squares add up to the least. That is least squares, as in Grade 12’s Misses and squares. NumPy’s np.polyfit(x, y, 1) does it and gives the line’s slope and intercept.

The points are each magnitude from 4.5 up, in steps of 0.1, with the logarithm of its count. A count under 10 is left out: when only three earthquakes are of 7.0 or more, one more or one fewer, by chance, moves the point a long way.

import numpy as npimport pandas as pd quakes = pd.read_csv("nepal_quakes.csv")quakes["year"] = quakes["date"].str[:4].astype(int)recent = quakes[quakes["year"] >= 1973]years = 53 + 10 / 12                 # January 1973 to the end of October 2026 mags = [round(4.5 + 0.1 * i, 1) for i in range(34)]          # 4.5, 4.6, … 7.8counts = [(recent["mag"] >= m).sum() for m in mags]use = [(m, n) for m, n in zip(mags, counts) if n >= 10]     # too few is chancex = [m for m, n in use]y = [np.log10(n) for m, n in use] slope, intercept = np.polyfit(x, y, 1)print(f"fitted from {x[0]} to {x[-1]}: log10(count) = {intercept:.2f} {slope:.2f} × magnitude")print(f"each step of one magnitude: {10 ** -slope:.1f} times fewer") for m in [6, 7, 8]:    per_year = 10 ** (intercept + slope * m) / years    print(f"{m} or more: one every {1 / per_year:.0f} years")

fitted from 4.5 to 6.3: log10(count) = 7.55 -1.05 × magnitudeeach step of one magnitude: 11.3 times fewer6 or more: one every 3 years7 or more: one every 37 years8 or more: one every 415 years

The slope is −1.05: the b-value of Nepal and its edges is 1.05, and each step of one magnitude has about 11 times fewer earthquakes. From the line, 10 ** (intercept + slope * m) is how many earthquakes of m or more the line expects in those 53 years; divided by the years, how many a year.

Step 5: the whole world

Is the line only Nepal’s? world_quakes.csv counts every earthquake of each magnitude from 6.0 up, anywhere, since 1973. Adding up from the bottom of the table, the biggest first, gives “this big or bigger”: [::-1] turns a column upside down, cumsum() adds as it goes, and the second [::-1] turns it back.

import numpy as npimport pandas as pd world = pd.read_csv("world_quakes.csv")world["at_least"] = world["count"][::-1].cumsum()[::-1]   # add up from the biggestfor m in [6.0, 7.0, 8.0, 9.0]:    print(m, "or more:", world[world["mag"] == m]["at_least"].iloc[0]) use = world[world["mag"] <= 7.5]slope, intercept = np.polyfit(use["mag"], np.log10(use["at_least"]), 1)print(f"fitted from 6.0 to 7.5: each step of one magnitude, {10 ** -slope:.1f} times fewer")for m in [8.0, 9.0]:    print(f"{m} or more: the line says {10 ** (intercept + slope * m):.0f}")

6.0 or more: 74417.0 or more: 7418.0 or more: 409.0 or more: 2fitted from 6.0 to 7.5: each step of one magnitude, 10.0 times fewer8.0 or more: the line says 749.0 or more: the line says 7

From 6 to 7, the world’s count falls from 7,441 to 741: ten times, almost exactly. The line fitted from 6.0 to 7.5 says 10.0. Above that the real counts fall faster than the line: 40 earthquakes of magnitude 8 or more where it says 74, two of 9 or more (Sumatra in 2004 and Japan in 2011) where it says 7.

Your program found the same rule in two very different piles of data: in Nepal’s 1,511 earthquakes since 1973, about 11 times fewer for each step of one magnitude; in the whole world’s, about 10 times fewer. Nobody planned it, and no single fault makes it: it is what millions of breaks in the rock, of every size, add up to. That is why the line is useful. The small earthquakes are many, and a few decades of them say how often the rare big ones should come: in this part of the world, by your line, one of magnitude 6 or more about every 3 years, 7 or more about every 37 years, 8 or more about every 400 years.

And here is what the line can’t do.

• It can’t say when. “One every 37 years” is an average over a long time, not a timetable. Since 1973 there have been three earthquakes of 7 or more here, not the line’s one or two: two of them in 2015, 17 days apart, because a big earthquake sets off others near it. The rate also doesn’t grow with the years since the last one: an earthquake isn’t “due” as a bus is.
• It can’t be trusted far beyond the data. Your line was fitted from 4.5 to 6.3, and asked about magnitude 8. In the world’s data the line overestimated the biggest earthquakes; in the Himalaya, the record of great earthquakes over the last 800 years (in Can We Tell When Disaster Will Strike?) is the better guide to the biggest.
• It can’t say how much harm. The same magnitude does very different harm in a city of brick houses and in a forest, at 10 km deep or at 100 km.

Three things to try

• Earthquakes come in families. In Step 3, count only 2015, then every year but 2015. How much of the decade’s count was Gorkha’s family? Seismologists take out aftershocks before fitting the line, which is the ★★★ challenge below.
• Choose the start. Change the smallest magnitude in Step 4 from 4.5 to 4.0, then to 5.0. How much does the slope move? Change the first year to 1990, or to 2000.
• A second way to find b. The average magnitude of the earthquakes above 4.5 says how steep the line is: b is about 0.434 ÷ (the average − 4.45). (0.434 is log₁₀ of the number e, and 4.45 is 4.5 less half of a 0.1 step.) Work it out with recent[recent["mag"] >= 4.5]["mag"].mean(), and compare it with your slope. This way uses every earthquake equally, while the line fit gives each magnitude one point.

An open question

Big earthquakes are rare, so the line is fitted mostly on small ones. Does the same rule hold all the way up, on every fault, or does a fault such as the Himalaya’s have a biggest earthquake of its own, which the line can’t see? Nobody knows for certain: there are not yet enough great earthquakes in the record to say.

Three challenges, each checked on the page: write your function in the program above the checker (press Edit), then press Check my program. The examples show what your function should give; there are more tests you can’t see. In the first two, mags is a plain Python list of magnitudes, like the mag column.

★ How many this big or bigger?

Write count_at_least(mags, smallest): how many of the magnitudes in the list are smallest or more.

def count_at_least(mags, smallest):    # How many magnitudes in the list are smallest or more?    return 0 print(count_at_least([4.2, 5.0, 6.1, 4.9, 5.0], 5.0))

★★ Fit the line yourself

Write slope(xs, ys): the slope of the least-squares line through the points, rounded to 2 decimals, without NumPy. The slope is the sum of (x − x̄)(y − ȳ) over the sum of (x − x̄)², where x̄ and ȳ are the averages of the xs and of the ys.

def slope(xs, ys):    # The slope of the least-squares line, rounded to 2 decimals.    return 0 print(slope([1, 2, 3, 4], [3, 5, 7, 9]))

★★★ Take out the aftershocks

Each earthquake is [day, mag, lat, lon], in time order. Call an earthquake an aftershock if some earthquake earlier in the list, at least as big, came 30 days or less before it (the same day counts) and 100 km or less from it. Write mainshocks(quakes): how many earthquakes in the list are not aftershocks. km gives the distance between two places.

import math def km(lat1, lon1, lat2, lon2):    # distance in km, near enough for places a few hundred km apart    dx = (lon2 - lon1) * 111.2 * math.cos(math.radians((lat1 + lat2) / 2))    dy = (lat2 - lat1) * 111.2    return math.hypot(dx, dy) def mainshocks(quakes):    # How many earthquakes in the list are not aftershocks?    return 0 print(mainshocks([[0, 6.0, 28.0, 84.0], [3, 4.5, 28.2, 84.3], [50, 4.8, 28.1, 84.1]]))

• Our World: Earthquakes and Plates: every earthquake of magnitude 6 or more since 1900 on a globe, played year by year, and the plates they mark.
• Can We Tell When Disaster Will Strike?: Nepal’s great earthquakes over 800 years, and the difference between a forecast and a prediction.
• Can You Design an Earthquake-Safe House?: what to do about the line’s answer.
• Grade 10, Data Science: pandas, from reading a table to groupby.
• Grade 12, Fitting a line: least squares, how well a line fits, and how sure its slope is.
• Which Way Would Rain Flow from Your Roof?: more real data from Nepal, its heights measured from space.