Formation Pressure Evaluation - Pore Pressure Prediction, Overpressure Mechanisms, and Drilling Mud Weight Window Design
Formation pressure evaluation is the discipline that determines where fluids in subsurface rocks are under pressure different from the normal hydrostatic gradient, and quantifies that pressure accurately enough to design the drilling mud weight that will control the well without fracturing the formation. The stakes are binary and catastrophic at both extremes: a mud weight too low allows formation fluids to influx into the wellbore, potentially leading to a blowout; a mud weight too high fractures the formation, causing lost circulation that can cost $500,000 to $5,000,000 per incident in offshore drilling and can destabilize the wellbore sufficiently to trap the drill string. The pressure window between these two extremes - the difference between pore pressure gradient and fracture gradient - can be as wide as 5 ppg in normally pressured onshore wells and as narrow as 0.3-0.5 ppg in deepwater wells with overpressured formations, where drilling a single section requires real-time pressure monitoring and mud weight adjustments measured in tenths of a ppg. The geomechanical understanding that enables accurate pore pressure prediction - from the mechanisms that create overpressure to the seismic velocity and resistivity signatures that reveal it before drilling - is one of the most commercially valuable technical skills in exploration and drilling engineering. This guide covers the complete pore pressure evaluation workflow: the physical mechanisms that generate overpressure, the quantitative prediction methods used pre-drill, the real-time monitoring techniques used while drilling, and the fracture gradient prediction methods that define the upper boundary of the mud weight window.
1. Overpressure Mechanisms - Why Formation Pressure Exceeds Hydrostatic
1.1 Normal Compaction and Hydrostatic Pressure
At any subsurface depth, the normal pore pressure is the hydrostatic pressure exerted by a continuous column of formation water from surface to that depth. When rocks compact normally during burial, pore fluids drain freely and pore pressure remains hydrostatic. Overpressure develops when some mechanism prevents this normal drainage, trapping fluids and building pressure above hydrostatic:
Normal hydrostatic pressure gradient calculation:
Pp_normal (psi) = rho_water x g x TVD = 0.433 x SG_water x TVD
Where SG_water = specific gravity of formation water (typically 1.02-1.10 for saline brine)
For SG = 1.07 (typical formation brine):
Pp_gradient = 0.433 x 1.07 = 0.463 psi/ft (8.9 ppg equivalent)
At TVD = 12,500 ft:
Pp_normal = 0.463 x 12,500 = 5,788 psi
Overburden (total vertical stress) gradient:
Sv (psi) = integral(rho_bulk x g dz) from surface to TVD
Average bulk density of overburden: 2.35 g/cm3 (for Gulf of Mexico deepwater)
Sv_gradient = 2.35 x 0.433 = 1.017 psi/ft → Sv at 12,500 ft = 12,713 psi
Effective stress (sigma_eff) = Sv - Pp:
Normal: sigma_eff = 12,713 - 5,788 = 6,925 psi (effective stress supports overburden)
Terzaghi effective stress principle:
When pore pressure increases toward Sv: sigma_eff → 0 → rock loses mechanical strength
At Pp = Sv (fracture gradient): sigma_eff = 0 → formation cannot support overburden → will fracture
This is why fracture gradient ≤ Sv everywhere in the subsurface
1.2 The Five Principal Overpressure Mechanisms
| Mechanism | Physical Process | Seismic/Log Signature | Pressure Magnitude |
|---|---|---|---|
| Disequilibrium compaction (undercompaction) | Rapid burial of low-permeability shales prevents pore water drainage. Rock compacts without dewatering → pore pressure builds above hydrostatic. Most common mechanism globally. | High sonic DT (slow velocity) relative to normal compaction trend. High neutron porosity. Low density. Seismic velocity decreases below expected trend. | Moderate to high: 0.6-0.85 psi/ft |
| Hydrocarbon generation and cracking | Kerogen-to-oil conversion generates fluids of larger molar volume. Oil-to-gas cracking at high maturity generates large volume increases (1 volume oil → 3-5 volumes gas). Volumetric expansion generates pressure in sealed system. | Anomalous resistivity increase (gas effect). Amplitude anomalies on seismic. Coincides with source rock maturation zones. | Can reach near-lithostatic: 0.85-0.95 psi/ft |
| Fluid expansion (aquathermal pressuring) | Pore fluids expand with increasing temperature during burial. If system is sealed, thermal expansion generates pressure. Significant only in very low-permeability sealed systems. | Difficult to distinguish from other mechanisms. Typically a secondary contributor. | Minor contribution: 0.1-0.2 psi/ft above hydrostatic |
| Clay diagenesis (smectite-illite transformation) | Smectite converts to illite at 80-120°C, releasing interlayer water that was bound in clay structure. Released water volume can exceed available pore space → pressure increase in sealed system. | Coincides with 80-120°C isotherm. Illite abundance from XRD increases simultaneously. Resistivity may decrease as conductive illite replaces non-conductive smectite. | Moderate: 0.6-0.75 psi/ft |
| Tectonic loading and lateral stress transfer | Horizontal compressional tectonics increase lateral stress on pore fluids. Thrust faulting transfers stress from deeper to shallower levels. Common in fold-and-thrust belts. | Pressure appears anomalously high for burial depth. No corresponding velocity anomaly (rock may still be normally compacted). Drilling confirms pressure without seismic prediction. | Variable: 0.7-0.90 psi/ft |
2. Pore Pressure Prediction Methods
2.1 Eaton's Method - The Industry Standard
Eaton's method (1975) is the most widely used quantitative pore pressure prediction technique. It estimates pore pressure from the deviation of sonic velocity (or resistivity) from the expected normal compaction trend (NCT). The physical basis is simple: under-compacted shales retain high porosity and slow velocity because their pore pressure prevented normal mechanical compaction. The deviation from the NCT quantifies the degree of under-compaction and therefore the excess pore pressure:
Eaton's method using sonic log:
Pp = Sv - (Sv - Pp_normal) x (DT_normal/DT_observed)^3
Where:
Pp = predicted pore pressure (psi or psi/ft gradient)
Sv = overburden stress (psi)
Pp_normal = hydrostatic pore pressure (psi)
DT_normal = sonic DT on the Normal Compaction Trend (NCT) at that depth (μs/ft)
DT_observed = measured sonic DT at that depth (μs/ft)
Exponent = 3 (Eaton's original, sometimes adjusted for local calibration)
Normal Compaction Trend (NCT) establishment:
Plot DT vs depth for shale intervals only (filter out sands, carbonates, volcanics)
In normally pressured intervals, DT decreases linearly with depth on a semi-log plot:
DT_normal(z) = DT_surface x exp(-k x z)
Where DT_surface ≈ 200 μs/ft and k = compaction coefficient (typically 0.0001-0.0002 ft^-1)
For Gulf of Mexico deepwater: DT_normal(z) = 200 x exp(-0.000125 x z)
At z = 12,500 ft: DT_normal = 200 x exp(-0.000125 x 12,500) = 200 x exp(-1.5625) = 200 x 0.2096 = 41.9 μs/ft
Eaton pore pressure calculation at 12,500 ft TVD:
Sv = 12,713 psi (from overburden integration above)
Pp_normal = 5,788 psi
DT_normal = 41.9 μs/ft
DT_observed (from sonic log in overpressured shale) = 68 μs/ft
Pp = 12,713 - (12,713 - 5,788) x (41.9/68)^3
= 12,713 - 6,925 x (0.6162)^3
= 12,713 - 6,925 x 0.2341
= 12,713 - 1,621 = 11,092 psi predicted pore pressure
Pore pressure gradient: 11,092/12,500 = 0.887 psi/ft = 17.1 ppg equivalent mud weight (EMW)
Normal: 0.463 psi/ft (8.9 ppg)
Predicted: 0.887 psi/ft (17.1 ppg)
Overpressure = 0.887 - 0.463 = 0.424 psi/ft above normal (8.2 ppg above normal hydrostatic)
This 17.1 ppg prediction is the pre-drill estimate. Actual mud weight must be 17.1 + 0.3 ppg safety margin = 17.4 ppg minimum mud weight to maintain wellbore control.
2.2 Eaton's Method Using Resistivity
Eaton resistivity method:
Pp = Sv - (Sv - Pp_normal) x (Ro/Rn)^1.2
Where:
Ro = measured shale resistivity at depth of interest (ohm·m)
Rn = normal compaction trend resistivity at same depth (ohm·m)
Exponent = 1.2 (Eaton's original for resistivity)
Normal compaction trend for resistivity:
Resistivity increases with depth in normally compacted shale (decreasing porosity → increasing resistivity)
Rn(z) = R_surface x exp(+k_R x z) (increasing exponential, opposite to velocity)
For Gulf of Mexico: Rn(z) = 1.2 x exp(0.000095 x z) ohm·m
At z = 12,500 ft: Rn = 1.2 x exp(0.000095 x 12,500) = 1.2 x exp(1.1875) = 1.2 x 3.279 = 3.935 ohm·m
Measured shale resistivity in overpressured zone: Ro = 1.85 ohm·m (lower than NCT → under-compacted → overpressured)
Pp_resistivity = 12,713 - (12,713 - 5,788) x (1.85/3.935)^1.2
= 12,713 - 6,925 x (0.4702)^1.2
0.4702^1.2: ln(0.4702) = -0.7551, -0.7551 x 1.2 = -0.9061, e^-0.9061 = 0.4042
= 12,713 - 6,925 x 0.4042
= 12,713 - 2,799 = 9,914 psi
Gradient: 9,914/12,500 = 0.793 psi/ft = 15.3 ppg
Sonic predicts 17.1 ppg, Resistivity predicts 15.3 ppg → discrepancy of 1.8 ppg
Resolution: The methods have different sensitivities to different overpressure mechanisms:
Sonic velocity is sensitive to BOTH compaction-related (mechanical) overpressure AND fluid property changes
Resistivity is sensitive primarily to compaction-related overpressure
If sonic predicts higher pressure than resistivity: likely a secondary mechanism (hydrocarbon generation, aquathermal) is contributing additional pressure beyond mechanical under-compaction
Design mud weight: Use the higher prediction (sonic, 17.1 ppg) for safety → MW = 17.4 ppg minimum
3. Fracture Gradient Prediction - The Upper Mud Weight Boundary
3.1 Hubbert and Willis / Matthews and Kelly Method
The fracture gradient is the pressure at which the formation will accept injected fluid and develop hydraulic fractures. It defines the maximum mud weight that can be used without inducing lost circulation. Exceeding the fracture gradient is a serious event in deepwater drilling where the narrow mud weight window may leave less than 0.5 ppg of margin between the minimum mud weight (pore pressure + safety margin) and the maximum mud weight (fracture gradient):
Eaton's fracture gradient equation (most commonly used):
FG = (nu/(1-nu)) x (Sv - Pp)/TVD + Pp/TVD
Where:
FG = fracture gradient (psi/ft or ppg)
nu = Poisson's ratio of the formation (dimensionless, 0.25-0.45 for shale)
Sv = overburden stress (psi)
Pp = pore pressure (psi)
TVD = true vertical depth (ft)
Fracture gradient calculation at 12,500 ft with overpressure:
nu = 0.42 (soft, young Gulf of Mexico shale)
Sv = 12,713 psi
Pp = 11,092 psi (Eaton sonic prediction)
FG = (0.42/(1-0.42)) x (12,713 - 11,092)/12,500 + 11,092/12,500
= (0.42/0.58) x 1,621/12,500 + 0.8874
= 0.7241 x 0.1297 + 0.8874
= 0.09394 + 0.8874 = 0.9813 psi/ft = 18.9 ppg fracture gradient
Available mud weight window:
Minimum MW (pore pressure + 0.3 ppg safety): 17.1 + 0.3 = 17.4 ppg
Maximum MW (fracture gradient - 0.5 ppg safety): 18.9 - 0.5 = 18.4 ppg
Available MW window = 18.4 - 17.4 = 1.0 ppg
At 12,500 ft in this well, the entire mud weight window is only 1.0 ppg wide. This is an extremely narrow window that requires:
1. Precise cuttings transport analysis (excess cuttings loading increases ECD above MW)
2. Managed Pressure Drilling (MPD) to control wellbore pressure in real-time
3. Continuous monitoring of flow-in vs flow-out for early influx/loss detection
Equivalent Circulating Density (ECD) check:
ECD = MW + annular pressure losses / (0.052 x TVD)
Annular pressure loss at 600 gpm in 8.5" hole: approximately 350 psi
ECD = 17.4 + 350/(0.052 x 12,500) = 17.4 + 0.538 = 17.94 ppg ECD while circulating
17.94 ppg ECD < 18.4 ppg maximum → Acceptable margin = 0.46 ppg while circulating
If reducing circulation rate to 400 gpm: annular loss ≈ 165 psi → ECD = 17.4 + 0.254 = 17.65 ppg → 0.75 ppg margin
3.2 LOT and XLOT Data - Direct Fracture Gradient Measurement
Leak-Off Test (LOT) procedure and interpretation:
The LOT is conducted after drilling out a new casing shoe to directly measure the fracture initiation pressure of the formation immediately below the shoe. It is the most reliable calibration point for fracture gradient prediction in subsequent sections:
LOT execution:
1. Drill 10-20 ft below casing shoe
2. Pull back to shoe and circulate to static conditions
3. Close annular BOP and inject mud at low rate (typically 0.25-0.5 bbl/min)
4. Record surface pressure vs cumulative volume pumped
5. Plot pressure vs volume → identify Leak-Off Point (LOP)
Extended LOT (XLOT) adds:
6. Continue pumping past LOP to Fracture Propagation Pressure (FPP)
7. Shut in and record Instantaneous Shut-In Pressure (ISIP) → Sh_min
8. Bleed off and re-pressurize → Fracture Reopening Pressure (FRP) → confirms Sh_min
Example LOT data at 9,500 ft casing shoe (6,000 ft TVD, 3,500 ft water depth):
Water depth converts: 3,500 ft water = 3,500 x 0.433 x 1.025 = 1,554 psi water pressure
TVD below mudline: 6,000 - 3,500 = 2,500 ft formation
Current mud weight: 10.8 ppg
Hydrostatic at shoe: 10.8 x 0.052 x 6,000 = 3,370 psi bottomhole pressure (BHP)
LOT conducted: Pumped 8 bbls before pressure deviation from linear
LOP surface pressure reading: 480 psi
BHP at LOP = 3,370 + 480 = 3,850 psi**
**Correction: in deepwater, measure from rotary kelly bushing (RKB), not from sea surface. RKB is 50 ft above sea surface → add 50 ft air gap correction:
Hydrostatic at shoe = 10.8 x 0.052 x (6,000 + 50) - 1.0 x 0.052 x 50 = 3,402 - 2.6 = 3,399 psi
BHP at LOP = 3,399 + 480 = 3,879 psi** (corrected)
Fracture gradient from LOT:
FG = BHP_LOP / TVD = 3,879 / 6,000 = 0.647 psi/ft = 12.44 ppg (at casing shoe depth)
ISIP (from XLOT): 380 psi surface → BHP = 3,399 + 380 = 3,779 psi
Sh_min = 3,779/6,000 = 0.630 psi/ft = 12.1 ppg minimum horizontal stress
This Sh_min directly calibrates the Eaton fracture gradient prediction for the next drill section below this shoe.
4. Real-Time Pore Pressure Monitoring While Drilling
4.1 MWD/LWD Pressure Monitoring and Drilling Exponent
Pre-drill pore pressure predictions are based on seismic velocity and offset well data. While drilling, real-time measurements update the pressure prediction and provide early warning of approaching overpressured zones. The drilling exponent (dc-exponent) is the oldest and still widely used real-time indicator of formation pressure, based on the observation that penetration rate changes systematically in overpressured formations:
Corrected d-exponent (dc-exponent) calculation:
dc = [log(ROP/60/N)] / [log(12 x WOB/(10^6 x D))] x (Rho_normal/Rho_mud)
Where:
ROP = Rate of penetration (ft/hr)
N = Rotary speed (RPM)
WOB = Weight on bit (lbs)
D = Bit diameter (inches)
Rho_normal = normal formation water gradient (ppg)
Rho_mud = actual mud weight (ppg)
dc-exponent calculation at two depths:
Depth 11,000 ft (normally pressured zone):
ROP=45 ft/hr, N=120 RPM, WOB=45,000 lbs, D=8.5", MW=12.5 ppg, Rho_normal=8.9 ppg
dc = [log(45/60/120)] / [log(12 x 45,000/(10^6 x 8.5))] x (8.9/12.5)
= [log(0.00625)] / [log(0.06353)] x 0.712
= [-2.2041] / [-1.1968] x 0.712
= 1.8417 x 0.712 = 1.311 (normal compaction trend value at 11,000 ft)
Depth 12,000 ft (approaching overpressured zone):
ROP=82 ft/hr (increased!), N=120 RPM, WOB=45,000 lbs, D=8.5", MW=12.5 ppg
dc = [log(82/60/120)] / [log(12 x 45,000/(10^6 x 8.5))] x (8.9/12.5)
= [log(0.01139)] / [-1.1968] x 0.712
= [-1.9434] / [-1.1968] x 0.712
= 1.6238 x 0.712 = 1.156
dc trend interpretation:
11,000 ft: dc = 1.311 (on NCT → normal pressure)
12,000 ft: dc = 1.156 (below NCT by 0.155 → formation is softer than expected → OVERPRESSURE WARNING)
The dc-exponent has DECREASED from 1.311 to 1.156 while the NCT predicts it should have INCREASED slightly with depth. This reversal in trend signals that the formation at 12,000 ft is under-compacted relative to normal → pore pressure is higher than expected.
Eaton dc-exponent pore pressure prediction:
Pp = Sv - (Sv - Pp_normal) x (dc_observed/dc_normal)^1.2
dc_normal at 12,000 ft: from NCT extrapolation = 1.380
Pp = 12,220 - (12,220 - 5,560) x (1.156/1.380)^1.2
= 12,220 - 6,660 x (0.8377)^1.2
0.8377^1.2: ln(0.8377)=-0.1772, -0.1772 x 1.2=-0.2126, e^-0.2126=0.8084
= 12,220 - 6,660 x 0.8084 = 12,220 - 5,384 = 6,836 psi → 0.570 psi/ft = 10.96 ppg
Current MW = 12.5 ppg → adequate safety margin maintained. But: trend is worsening.
Action: Plan contingent mud weight increase to 13.5 ppg within next 500 ft of drilling.
4.2 PWD (Pressure While Drilling) - Direct Annular Pressure Measurement
| Real-Time Indicator | What It Measures | Overpressure Warning Signal | Response Time |
|---|---|---|---|
| PWD annular pressure | Direct annular pressure at bit depth using MWD pressure sensor. Measures ECD continuously during drilling. | ECD decreasing without mud weight change → formation taking fluid (loss). ECD increasing without flow rate change → wellbore taking on gas (influx). | Real-time (seconds) |
| Flow check (pit gain) | Measure flow-out vs flow-in volume. Pit gain = more fluid returning than pumped = formation influx (kick). | Any pit gain > 0.5 bbls requires immediate well shut-in investigation. 2 bbls gain = definite kick in progress. | Minutes (continuous monitoring) |
| Gas readings (total gas, chromatograph) | Gas detector on mud returns measures C1-C5 hydrocarbons from cutting gases. Connection gas (spike on each connection) indicates near-balanced conditions. | Connection gas increasing in trend → mud weight approaching or below pore pressure. Background gas > 2x previous level → formation pressure increasing. | Minutes to hours (lag time) |
| Shale cuttings analysis | Cuttings morphology changes in overpressured shale: splintery/angular fragments, cavings from borehole wall sloughing, shale "popcorn" texture. | Large cavings arriving at surface → wellbore instability → either MW too low (plastic failure) or too high (tensile failure). Distinguish by caving morphology. | Hours (lag time to surface) |
5. Kick Detection and Well Control - The Emergency Response
5.1 Kick Volume and Pressure Calculations
Kick intensity and shut-in pressure analysis:
When a kick is detected (pit gain, flow with pumps off), the well is shut in and the following pressures are recorded:
SIDPP = Shut-In Drill Pipe Pressure (psi) - reflects pore pressure of influxed formation
SICP = Shut-In Casing Pressure (psi) - reflects SIDPP + gas migration effect
Formation pore pressure from SIDPP:
Pp = SIDPP + (MW x 0.052 x TVD)
Example: Kick detected at TVD = 12,500 ft
Current MW = 17.4 ppg (just above Eaton prediction of 17.1 ppg)
After shut-in: SIDPP = 285 psi
Pp = 285 + (17.4 x 0.052 x 12,500) = 285 + 11,310 = 11,595 psi actual pore pressure
Pp gradient = 11,595/12,500 = 0.928 psi/ft = 17.84 ppg actual EMW
Eaton predicted 17.1 ppg → actual is 17.84 ppg → prediction was low by 0.74 ppg
Kill mud weight required = 17.84 + 0.3 = 18.14 ppg (within the 18.4 ppg fracture gradient limit)
Kick volume and influx type determination:
Pit gain = 8 bbls
SIDPP = 285 psi, SICP = 485 psi
Difference SICP - SIDPP = 200 psi
Influx density calculation:
rho_influx = MW - (SICP - SIDPP) / (0.052 x kick_length_in_annulus)
Kick length in annulus (8.5" hole, 5" DP):
Annular capacity = (8.5^2 - 5^2)/1029.4 = (72.25-25)/1029.4 = 47.25/1029.4 = 0.04590 bbl/ft
Kick length = 8 bbls / 0.04590 = 174 ft
rho_influx = 17.4 - 200/(0.052 x 174) = 17.4 - 200/9.048 = 17.4 - 22.10 = -4.7 ppg
Negative density → clearly gas kick (gas has near-zero density, hence very large negative result)
Confirm: Gas kick identified → requires careful kill procedure with gas migration monitoring
Driller's method kill mud volume:
Kill MW = 18.14 ppg (from above)
Drill pipe volume to circulate kill mud to bit = dp_capacity x bit_depth
DP capacity = 5" x 4.276" ID = (4.276)^2/1029.4 = 0.01778 bbl/ft
Volume = 0.01778 x 12,500 = 222.3 bbls kill mud to bit
Conclusion
The Eaton sonic versus resistivity pore pressure discrepancy in this article - 17.1 ppg from sonic versus 15.3 ppg from resistivity at the same depth - is not a calculation error but a physically meaningful signal. When sonic velocity predicts higher pore pressure than resistivity for the same depth in the same shale, it is diagnostic evidence that a secondary overpressure mechanism (hydrocarbon generation or fluid expansion) is contributing pressure beyond the mechanical under-compaction that both methods equally detect. Resistivity responds primarily to porosity changes (compaction state), while sonic velocity responds to both porosity changes and the elastic moduli of the pore fluid (gas dramatically slows velocity, water does not). The divergence between the two predictions is therefore a direct indicator of the presence of gas or light hydrocarbons in the overpressured formation. Using the higher prediction for mud weight design is correct - but the divergence itself also tells the drilling team that they are approaching an interval with possible gas-charged overpressured shales, which changes the well control risk profile and the contingency planning requirements.
The kick analysis calculation - SIDPP of 285 psi identifying an actual pore pressure of 17.84 ppg against a pre-drill prediction of 17.1 ppg, with the influx density calculation confirming gas rather than water or oil - demonstrates the complete well control decision chain in action. The 0.74 ppg prediction error is not catastrophic: the actual pore pressure (17.84 ppg) is still below the fracture gradient (18.9 ppg), and the required kill mud weight (18.14 ppg) fits within the available mud weight window (18.4 ppg maximum). But the margin is now 0.26 ppg between the kill mud weight and the fracture gradient - a margin that evaporates quickly if there are cuttings surcharges or swab/surge effects during the kill operation. The drilling engineer managing this well control event must simultaneously track SICP for gas migration (gas migrates up at 500-1,000 ft/hr if well is shut in), plan pump rates that keep ECD below fracture gradient, and prepare for the possibility that the kill mud weight itself will induce lost circulation before the kick is fully killed.
For drilling engineers and geoscientists building expertise in pore pressure prediction and well control, the following references provide the essential framework: Pore Pressure Prediction and Drilling Geomechanics covers Eaton's method, Normal Compaction Trend establishment, and real-time pressure monitoring techniques, while Well Control and Blowout Prevention provides the quantitative well control calculations for kick detection, influx analysis, and kill procedures.
Want to access our pore pressure toolkit with Eaton sonic and resistivity calculator, NCT establishment tool, fracture gradient predictor, ECD calculator, kick intensity analyzer, and kill mud weight determination, or discuss pore pressure prediction for a specific well? Join our Telegram group for drilling engineering and pore pressure discussions, or visit our YouTube channel for step-by-step tutorials on Eaton's method, LOT interpretation, and real-time pressure monitoring.
Disclosure: This article contains affiliate links. If you purchase a book through these links, Petrosmart Academy may earn a small commission at no additional cost to you. We only recommend resources that provide genuine value to drilling engineers and geoscientists.

0 Comments