0
# Gravimetry and Subsidence
1

2
Gravity and subsidence modeling capabilities for time-lapse reservoir monitoring and geomechanical analysis, enabling integrated reservoir characterization through geophysical methods.
3

4
## Capabilities
5

6
### Gravity Modeling
7

8
Gravity simulation operations for time-lapse gravimetric monitoring of reservoir changes.
9

10
```python { .api }
11
class ResdataGrav:
12
    """Gravity simulation operations for reservoir monitoring."""
13
    
14
    def __init__(self, grid: Grid, porv_kw: ResdataKW):
15
        """
16
        Initialize gravity calculator.
17
        
18
        Args:
19
            grid (Grid): Reservoir grid
20
            porv_kw (ResdataKW): Pore volume keyword
21
        """
22
    
23
    def add_survey_GRAV(self, survey_name: str, x: float, y: float, z: float):
24
        """
25
        Add gravity survey point.
26
        
27
        Args:
28
            survey_name (str): Name of the survey
29
            x, y, z (float): Survey point coordinates
30
        """
31
    
32
    def add_survey_GRAVDIFF(self, survey_name: str, base_survey: str, 
33
                           x: float, y: float, z: float):
34
        """
35
        Add differential gravity survey point.
36
        
37
        Args:
38
            survey_name (str): Name of the survey
39
            base_survey (str): Reference survey name
40
            x, y, z (float): Survey point coordinates
41
        """
42
    
43
    def eval_grav(self, time: datetime = None) -> numpy.ndarray:
44
        """
45
        Evaluate gravity response.
46
        
47
        Args:
48
            time (datetime, optional): Evaluation time
49
            
50
        Returns:
51
            numpy.ndarray: Gravity values in milligals (mGal)
52
        """
53
    
54
    def load_PORV(self, porv_kw: ResdataKW):
55
        """Load pore volume data."""
56
    
57
    def load_RHO(self, rho_kw: ResdataKW):
58
        """Load density data."""
59
    
60
    def add_std_density(self, oil_density: float, gas_density: float, 
61
                       water_density: float):
62
        """
63
        Add standard fluid densities.
64
        
65
        Args:
66
            oil_density (float): Oil density (kg/m³)
67
            gas_density (float): Gas density (kg/m³)
68
            water_density (float): Water density (kg/m³)
69
        """
70
```
71

72
### Subsidence Modeling
73

74
Subsidence simulation operations for geomechanical analysis and surface deformation monitoring.
75

76
```python { .api }
77
class ResdataSubsidence:
78
    """Subsidence simulation operations for geomechanical analysis."""
79
    
80
    def __init__(self, grid: Grid):
81
        """
82
        Initialize subsidence calculator.
83
        
84
        Args:
85
            grid (Grid): Reservoir grid
86
        """
87
    
88
    def add_survey_SUBSIDENCE(self, survey_name: str, x: float, y: float, z: float):
89
        """
90
        Add subsidence monitoring point.
91
        
92
        Args:
93
            survey_name (str): Name of the survey
94
            x, y, z (float): Monitoring point coordinates
95
        """
96
    
97
    def eval_subsidence(self, time: datetime = None) -> numpy.ndarray:
98
        """
99
        Evaluate subsidence response.
100
        
101
        Args:
102
            time (datetime, optional): Evaluation time
103
            
104
        Returns:
105
            numpy.ndarray: Subsidence values in meters
106
        """
107
    
108
    def load_PRESSURE(self, pressure_kw: ResdataKW):
109
        """Load pressure data for geomechanical calculations."""
110
    
111
    def load_PORV(self, porv_kw: ResdataKW):
112
        """Load pore volume data."""
113
    
114
    def add_node(self, x: float, y: float, z: float):
115
        """
116
        Add calculation node.
117
        
118
        Args:
119
            x, y, z (float): Node coordinates
120
        """
121
```
122

123
## Usage Examples
124

125
### Gravity Survey Analysis
126

127
```python
128
from resdata.gravimetry import ResdataGrav
129
from resdata.grid import Grid
130
from resdata.resfile import ResdataFile
131
import numpy as np
132
import matplotlib.pyplot as plt
133

134
# Load grid and simulation data
135
grid = Grid("SIMULATION.EGRID")
136
init_file = ResdataFile("SIMULATION.INIT")
137
restart_file = ResdataFile("SIMULATION.UNRST")
138

139
# Get pore volume
140
porv = init_file.get_kw("PORV")
141

142
# Initialize gravity calculator
143
grav_calc = ResdataGrav(grid, porv)
144

145
# Set standard fluid densities (kg/m³)
146
grav_calc.add_std_density(
147
    oil_density=850.0,    # Light oil
148
    gas_density=200.0,    # Gas at reservoir conditions
149
    water_density=1020.0  # Formation water
150
)
151

152
# Define gravity survey grid
153
survey_points = []
154
x_range = np.linspace(0, 4000, 21)  # 21 points over 4 km
155
y_range = np.linspace(0, 3000, 16)  # 16 points over 3 km
156
z_survey = 0.0  # Surface elevation
157

158
for i, y in enumerate(y_range):
159
    for j, x in enumerate(x_range):
160
        survey_name = f"GRAV_{i:02d}_{j:02d}"
161
        grav_calc.add_survey_GRAV(survey_name, x, y, z_survey)
162
        survey_points.append((x, y, survey_name))
163

164
print(f"Added {len(survey_points)} gravity survey points")
165

166
# Evaluate gravity response
167
gravity_response = grav_calc.eval_grav()
168
print(f"Gravity response shape: {gravity_response.shape}")
169
print(f"Gravity range: {gravity_response.min():.3f} to {gravity_response.max():.3f} mGal")
170

171
# Reshape for plotting
172
nx, ny = len(x_range), len(y_range)
173
gravity_grid = gravity_response.reshape(ny, nx)
174

175
# Plot gravity map
176
plt.figure(figsize=(12, 9))
177
contour = plt.contourf(x_range, y_range, gravity_grid, levels=20, cmap='RdBu_r')
178
plt.colorbar(contour, label='Gravity Anomaly (mGal)')
179
plt.xlabel('X (m)')
180
plt.ylabel('Y (m)')
181
plt.title('Gravity Survey - Base Case')
182
plt.axis('equal')
183
plt.grid(True, alpha=0.3)
184
plt.show()
185
```
186

187
### Time-Lapse Gravity Analysis
188

189
```python
190
from resdata.gravimetry import ResdataGrav
191
from resdata.grid import Grid
192
from resdata.resfile import ResdataFile
193
import numpy as np
194
import matplotlib.pyplot as plt
195

196
# Load data
197
grid = Grid("SIMULATION.EGRID")
198
init_file = ResdataFile("SIMULATION.INIT")
199
restart_file = ResdataFile("SIMULATION.UNRST")
200

201
# Get base data
202
porv = init_file.get_kw("PORV")
203

204
# Initialize gravity calculator
205
grav_calc = ResdataGrav(grid, porv)
206
grav_calc.add_std_density(850.0, 200.0, 1020.0)
207

208
# Add survey points over reservoir
209
survey_points = []
210
for i in range(5):
211
    for j in range(5):
212
        x = 1000 + j * 500  # 500m spacing
213
        y = 1000 + i * 500
214
        z = 0.0
215
        survey_name = f"GRAV_{i}_{j}"
216
        grav_calc.add_survey_GRAV(survey_name, x, y, z)
217
        survey_points.append((x, y, survey_name))
218

219
# Evaluate gravity at different time steps
220
time_steps = [0, 365, 730, 1095]  # Days: initial, 1yr, 2yr, 3yr
221
gravity_maps = {}
222

223
for day in time_steps:
224
    # Load density data for this time step
225
    # (This would typically involve calculating densities from saturation data)
226
    
227
    # For demonstration, evaluate gravity response
228
    gravity_response = grav_calc.eval_grav()
229
    gravity_maps[day] = gravity_response.reshape(5, 5)
230
    
231
    print(f"Day {day}: Gravity range {gravity_response.min():.3f} to "
232
          f"{gravity_response.max():.3f} mGal")
233

234
# Plot time-lapse gravity changes
235
fig, axes = plt.subplots(2, 2, figsize=(15, 12))
236
axes = axes.flatten()
237

238
x_coords = [1000 + j * 500 for j in range(5)]
239
y_coords = [1000 + i * 500 for i in range(5)]
240

241
for idx, day in enumerate(time_steps):
242
    ax = axes[idx]
243
    im = ax.contourf(x_coords, y_coords, gravity_maps[day], 
244
                     levels=20, cmap='RdBu_r')
245
    ax.set_title(f'Day {day} ({day/365:.1f} years)')
246
    ax.set_xlabel('X (m)')
247
    ax.set_ylabel('Y (m)')
248
    ax.axis('equal')
249
    fig.colorbar(im, ax=ax, label='Gravity (mGal)')
250

251
plt.tight_layout()
252
plt.show()
253

254
# Calculate gravity changes
255
base_gravity = gravity_maps[0]
256
for day in time_steps[1:]:
257
    gravity_change = gravity_maps[day] - base_gravity
258
    max_change = np.max(np.abs(gravity_change))
259
    print(f"Maximum gravity change at day {day}: {max_change:.3f} mGal")
260
```
261

262
### Subsidence Monitoring
263

264
```python
265
from resdata.gravimetry import ResdataSubsidence
266
from resdata.grid import Grid
267
from resdata.resfile import ResdataFile
268
import numpy as np
269
import matplotlib.pyplot as plt
270

271
# Load grid and data
272
grid = Grid("SIMULATION.EGRID")
273
restart_file = ResdataFile("SIMULATION.UNRST")
274

275
# Initialize subsidence calculator
276
subsidence_calc = ResdataSubsidence(grid)
277

278
# Add monitoring points at surface
279
monitoring_points = []
280
for i in range(11):
281
    for j in range(11):
282
        x = j * 400  # 400m spacing over 4km
283
        y = i * 300  # 300m spacing over 3km
284
        z = 0.0      # Surface level
285
        
286
        subsidence_calc.add_survey_SUBSIDENCE(f"SUB_{i}_{j}", x, y, z)
287
        monitoring_points.append((x, y))
288

289
print(f"Added {len(monitoring_points)} subsidence monitoring points")
290

291
# Load pressure data
292
pressure = restart_file.get_kw("PRESSURE", 0)  # Initial pressure
293
subsidence_calc.load_PRESSURE(pressure)
294

295
# Get pore volume if available
296
try:
297
    porv = restart_file.get_kw("PORV")
298
    subsidence_calc.load_PORV(porv)
299
except:
300
    print("PORV not found in restart file")
301

302
# Evaluate subsidence
303
subsidence_response = subsidence_calc.eval_subsidence()
304
print(f"Subsidence range: {subsidence_response.min():.4f} to "
305
      f"{subsidence_response.max():.4f} m")
306

307
# Reshape for plotting
308
subsidence_grid = subsidence_response.reshape(11, 11)
309

310
# Plot subsidence map
311
x_coords = [j * 400 for j in range(11)]
312
y_coords = [i * 300 for i in range(11)]
313

314
plt.figure(figsize=(12, 9))
315
contour = plt.contourf(x_coords, y_coords, subsidence_grid, 
316
                      levels=20, cmap='RdYlBu')
317
plt.colorbar(contour, label='Subsidence (m)')
318
plt.xlabel('X (m)')
319
plt.ylabel('Y (m)')
320
plt.title('Surface Subsidence Map')
321
plt.axis('equal')
322
plt.grid(True, alpha=0.3)
323

324
# Add contour lines
325
contour_lines = plt.contour(x_coords, y_coords, subsidence_grid, 
326
                           levels=10, colors='black', alpha=0.4, linewidths=0.5)
327
plt.clabel(contour_lines, inline=True, fontsize=8, fmt='%.3f')
328

329
plt.show()
330

331
# Identify maximum subsidence
332
max_subsidence_idx = np.unravel_index(np.argmax(subsidence_grid), subsidence_grid.shape)
333
max_x = x_coords[max_subsidence_idx[1]]
334
max_y = y_coords[max_subsidence_idx[0]]
335
max_value = subsidence_grid[max_subsidence_idx]
336

337
print(f"Maximum subsidence: {max_value:.4f} m at ({max_x}, {max_y})")
338
```
339

340
### Integrated Gravity-Subsidence Analysis
341

342
```python
343
from resdata.gravimetry import ResdataGrav, ResdataSubsidence
344
from resdata.grid import Grid
345
from resdata.resfile import ResdataFile
346
import numpy as np
347
import matplotlib.pyplot as plt
348

349
# Load data
350
grid = Grid("SIMULATION.EGRID")
351
init_file = ResdataFile("SIMULATION.INIT")
352
restart_file = ResdataFile("SIMULATION.UNRST")
353

354
# Get required data
355
porv = init_file.get_kw("PORV")
356
pressure = restart_file.get_kw("PRESSURE", 0)
357

358
# Initialize both calculators
359
grav_calc = ResdataGrav(grid, porv)
360
grav_calc.add_std_density(850.0, 200.0, 1020.0)
361

362
subsidence_calc = ResdataSubsidence(grid)
363
subsidence_calc.load_PRESSURE(pressure)
364
subsidence_calc.load_PORV(porv)
365

366
# Define common survey grid
367
survey_grid = []
368
for i in range(7):
369
    for j in range(7):
370
        x = 500 + j * 500  # 500m spacing
371
        y = 500 + i * 500
372
        z = 0.0
373
        
374
        # Add to both calculators
375
        grav_calc.add_survey_GRAV(f"POINT_{i}_{j}", x, y, z)
376
        subsidence_calc.add_survey_SUBSIDENCE(f"POINT_{i}_{j}", x, y, z)
377
        
378
        survey_grid.append((x, y))
379

380
# Evaluate both responses
381
gravity_response = grav_calc.eval_grav()
382
subsidence_response = subsidence_calc.eval_subsidence()
383

384
# Reshape for plotting
385
gravity_grid = gravity_response.reshape(7, 7)
386
subsidence_grid = subsidence_response.reshape(7, 7)
387

388
# Plot combined analysis
389
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18, 5))
390

391
x_coords = [500 + j * 500 for j in range(7)]
392
y_coords = [500 + i * 500 for i in range(7)]
393

394
# Gravity map
395
im1 = ax1.contourf(x_coords, y_coords, gravity_grid, levels=15, cmap='RdBu_r')
396
ax1.set_title('Gravity Anomaly')
397
ax1.set_xlabel('X (m)')
398
ax1.set_ylabel('Y (m)')
399
fig.colorbar(im1, ax=ax1, label='Gravity (mGal)')
400
ax1.axis('equal')
401

402
# Subsidence map
403
im2 = ax2.contourf(x_coords, y_coords, subsidence_grid, levels=15, cmap='RdYlBu')
404
ax2.set_title('Surface Subsidence')
405
ax2.set_xlabel('X (m)')
406
ax2.set_ylabel('Y (m)')
407
fig.colorbar(im2, ax=ax2, label='Subsidence (m)')
408
ax2.axis('equal')
409

410
# Correlation plot
411
ax3.scatter(gravity_response, subsidence_response * 1000, alpha=0.7)  # Convert to mm
412
ax3.set_xlabel('Gravity Anomaly (mGal)')
413
ax3.set_ylabel('Subsidence (mm)')
414
ax3.set_title('Gravity vs Subsidence Correlation')
415
ax3.grid(True, alpha=0.3)
416

417
# Calculate correlation coefficient
418
correlation = np.corrcoef(gravity_response, subsidence_response)[0, 1]
419
ax3.text(0.05, 0.95, f'Correlation: {correlation:.3f}', 
420
         transform=ax3.transAxes, verticalalignment='top',
421
         bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
422

423
plt.tight_layout()
424
plt.show()
425

426
print(f"Gravity-Subsidence Analysis Summary:")
427
print(f"  Gravity range: {gravity_response.min():.3f} to {gravity_response.max():.3f} mGal")
428
print(f"  Subsidence range: {subsidence_response.min()*1000:.2f} to {subsidence_response.max()*1000:.2f} mm")
429
print(f"  Correlation coefficient: {correlation:.3f}")
430
```
431

432
### Survey Design Optimization
433

434
```python
435
from resdata.gravimetry import ResdataGrav
436
from resdata.grid import Grid
437
from resdata.resfile import ResdataFile
438
import numpy as np
439
import matplotlib.pyplot as plt
440

441
# Load data
442
grid = Grid("SIMULATION.EGRID")
443
init_file = ResdataFile("SIMULATION.INIT")
444
porv = init_file.get_kw("PORV")
445

446
# Initialize gravity calculator
447
grav_calc = ResdataGrav(grid, porv)
448
grav_calc.add_std_density(850.0, 200.0, 1020.0)
449

450
# Test different survey designs
451
survey_designs = {
452
    'coarse': {'spacing': 1000, 'description': '1km spacing'},
453
    'medium': {'spacing': 500, 'description': '500m spacing'},
454
    'fine': {'spacing': 250, 'description': '250m spacing'}
455
}
456

457
results = {}
458

459
for design_name, design in survey_designs.items():
460
    spacing = design['spacing']
461
    
462
    # Create new gravity calculator for this design
463
    grav_test = ResdataGrav(grid, porv)
464
    grav_test.add_std_density(850.0, 200.0, 1020.0)
465
    
466
    # Add survey points
467
    survey_points = []
468
    x_range = np.arange(0, 4001, spacing)
469
    y_range = np.arange(0, 3001, spacing)
470
    
471
    for i, y in enumerate(y_range):
472
        for j, x in enumerate(x_range):
473
            grav_test.add_survey_GRAV(f"P_{i}_{j}", x, y, 0.0)
474
            survey_points.append((x, y))
475
    
476
    # Evaluate gravity
477
    gravity_response = grav_test.eval_grav()
478
    
479
    results[design_name] = {
480
        'points': len(survey_points),
481
        'spacing': spacing,
482
        'gravity_range': gravity_response.max() - gravity_response.min(),
483
        'gravity_std': np.std(gravity_response),
484
        'response': gravity_response
485
    }
486

487
# Compare survey designs
488
print("Survey Design Comparison:")
489
print(f"{'Design':<10} {'Points':<8} {'Spacing':<10} {'Range':<12} {'Std Dev':<10}")
490
print("-" * 60)
491

492
for name, result in results.items():
493
    print(f"{name:<10} {result['points']:<8} {result['spacing']:<10} "
494
          f"{result['gravity_range']:<12.3f} {result['gravity_std']:<10.3f}")
495

496
# Plot survey design comparison
497
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
498

499
for idx, (name, result) in enumerate(results.items()):
500
    ax = axes[idx]
501
    
502
    spacing = result['spacing']
503
    x_range = np.arange(0, 4001, spacing)
504
    y_range = np.arange(0, 3001, spacing)
505
    
506
    gravity_grid = result['response'].reshape(len(y_range), len(x_range))
507
    
508
    im = ax.contourf(x_range, y_range, gravity_grid, levels=15, cmap='RdBu_r')
509
    ax.set_title(f"{name.title()} Grid ({result['points']} points)")
510
    ax.set_xlabel('X (m)')
511
    ax.set_ylabel('Y (m)')
512
    fig.colorbar(im, ax=ax, label='Gravity (mGal)')
513
    ax.axis('equal')
514

515
plt.tight_layout()
516
plt.show()
517

518
# Cost-benefit analysis
519
print("\nCost-Benefit Analysis:")
520
base_cost_per_point = 1000  # Cost per gravity measurement
521
for name, result in results.items():
522
    total_cost = result['points'] * base_cost_per_point
523
    resolution = result['gravity_std']
524
    efficiency = resolution / (total_cost / 1000)  # Resolution per k$
525
    
526
    print(f"{name:<10}: Cost ${total_cost:,}, "
527
          f"Resolution {resolution:.3f} mGal, "
528
          f"Efficiency {efficiency:.6f} mGal/k$")
529
```
530

531
## Types
532

533
```python { .api }
534
# Gravity and subsidence response arrays
535
GravityResponse = numpy.ndarray      # Gravity values in milligals
536
SubsidenceResponse = numpy.ndarray   # Subsidence values in meters
537

538
# Survey point coordinates
539
SurveyPoint = tuple[float, float, float]  # (x, y, z)
540
SurveyGrid = list[SurveyPoint]
541

542
# Fluid density parameters
543
FluidDensities = dict[str, float]    # Keys: 'oil', 'gas', 'water' (kg/m³)
544

545
# Geophysical survey metadata
546
SurveyMetadata = dict[str, any]      # Survey configuration and parameters
547
```