Method for Simulation and Optimization of Underground Gas Storage Performance

. Abstract — Method for Simulation and Optimization of Underground Gas Storage Performance — Proper simulation and identiﬁcation of the ﬂow potential of a gas storage plant is only possible if the nonlinear limits related to cavern operation and an optimal strategy (set of decision rules) related to the plant operation mode are considered. An efﬁcient operation of a storage plant is a challenging task due to the complexity of cavern rock mechanical restrictions, as well as other technical restrictions imposed on different plant devices. The scope of this paper is to cater for these challenges by deﬁning a set of different operational strategies, i.e., sets of actions which constitute the predeﬁned high level mechanism that allows for an economic and efﬁcient plant operation. In this paper, one speciﬁc strategy example is described in detail. Speciﬁcally, we give simulation results and outline an optimization procedure designed to maximize the plant’s performance capabilities. The general plant model and the form of rock mechanical restrictions of cavern operation are reviewed, and the construction of the computational algorithm is analyzed.


INTRODUCTION
Gas storage facilities, based either on underground salt caverns or depleted gas fields, are the key element in an industrial gas grid.Operation of the cavern plant, which is used not only for seasonal gas storing, but also for balancing on the daily or weekly scale, is a complex task.At any time point an operator can undertake many different actions, therefore the proper computational tool is required.The aim of our article is to build a framework that comprises all elements involved in the plant operation; having the model of caverns and above ground devices (compressors, preheaters, cyclon and glycol separators, etc.) does not in itself give a straightforward answer about how to operate the plant efficiently.The presented framework is intended to be used in order to choose the most optimal way of plant operation and to simulate the future gas storage states.The operational needs of industrial gas storage are determined by gas market demands, however secure plant operation is a major factor to be considered while fulfilling these demands.Due to the complex form of several different requirements, one of the most important parameters of a gas storage facility, i.e., maximal hourly injection and withdrawal gas flow, is sometimes artificially undervalued.All factors limiting the accuracy of the forecast of injection and withdrawal flow potential result in a quite conservative plant operation.The lack of an accurate prognosis therefore obstructs an efficient operation of the plant within the safety margins.The factors contributing to the difficulty of the forecasting task include: -nonlinear characteristics of caverns [1-3] and above ground facilities (including compressors); -complex nonlinear cavern operation constraints [4][5][6][7]; -the possibility of violating geological cavern constraints which may lead to the temporary closure of the cavern [8][9][10]; -scale of the gas storage facility which can include several caverns and several compressors [11]; -long term time horizon (some time dependent cavern operation restrictions have a one year cycle); -quick cycling of operation mode [12].
The most complex form of rock-mechanical restrictions is connected with the operation of the cavern in the low-pressure region, the analysis of such a case is discussed in [13][14][15].From the perspective of integration of all elements into the overall framework, the challenging task is to include the time-dependent form of rock mechanics restrictions.
In this paper, the Storage Operation Expert (SOE) method is presented.It takes into consideration the thermodynamic model of gas caverns, model of above ground facilities, nonlinear constraints of rock mechanics imposed on the cavern operation and a predefined strategy for utilization of caverns.This method can be used to support plant operation in terms of technical and commercial concerns.Depending on simulation inputs, this method provides information about: -realistic hourly plant withdrawal and injection limits for the flexible time horizon.This information can be easily converted into nomination limits available to storage users (based on individual storage users curves).To calculate hourly plant limits one should define an infinite (for injection) or negative infinite (for withdrawal) gas flow demand (refer Sect.3.2.);-state of rock mechanics timers that are active when a single cavern pressure is in a restricted range related with rules and limits of rock mechanics.Based on this information an operator can see how many days a cavern was operating in the restricted range.Specifically, the operator is informed about how many days he has to safely get the cavern out of the restricted range and how many days he can still safely operate in restricted ranges; -optimal realization of aggregated business nominations for the plant.This optimal realization is a set of trajectories for each cavern in the plant.Each trajectory provides information about the gas flow associated with each cavern for every step of simulation (e.g.every hour of plant operation).Moreover, an operator can obtain the information about trajectories of wellhead pressure, temperature and rock mechanic timers.In order to calculate this information, one should set the gas flow demand to the accepted aggregated business nominations trajectory.Furthermore, our method can be used to perform ''WHAT-IF'' simulations which allow operators to test the plant for different gas demand scenarios, different plant operation strategies, different availability and configuration of Above Ground Devices (AGDs) and different availability and configuration of gas caverns.The presented solution performs simulation with the resolution of one hour or one day.The simulator calculates an optimal decision trajectory for each cavern in the storage according to the following equation: where k denotes the cavern number and t is the time (in units of simulation steps).This decision trajectory determines the time evolution of calculated storage and cavern states.The cavern state is defined by: -wellhead pressure, volume and temperature; -timers and other information related with the rock mechanics limits, i.e., how long a cavern can operate below p LL (lower limit of unrestricted operation) pressure level, or how long a cavern still has to wait above p R, recreation pressure level (see Sect. 2); -connected/disconnected status.The storage state is defined by: -hourly or daily volume increment; -mode (injection with compressors, withdrawal with compressors, free flow injection or free flow withdrawal); -availability of caverns and AGDs; -configuration of AGDs (e.g. two out of three compressors running).

STORAGE PLANT MODEL
The plant model used in the simulations (Fig. 1) consists of five main parts: 1. system pipe, which connects a set of AGDs with the gas grid.The pressure in the system pipe is the same as the pressure in the gas grid; 2. set of AGDs.This module represents all AGDs which take part in the withdrawal or injection action; 3. field pipe which connects caverns with the set of AGDs.It is assumed that: during withdrawal, pressure in the field pipe is dp wdt (for example dp wdt ¼ 2 bar) lower than the lowest pressure among caverns connected to the field pipe.During injection, the pressure in the field pipe is dp inj (for example dp inj ¼ 3 bar) larger than the highest pressure among caverns connected to the field pipe; 4. pressure control valves.Individual caverns may be disconnected or connected to the field pipe.When they are connected, a dedicated controller controls the gas flow to and from the caverns by acting on the control valves.Due to the control operation range and safety regulations, it is assumed that the difference between the maximum and minimum pressure in the set of connected caverns cannot exceed p window (for example p window ¼ 50 bar); 5. gas caverns are modeled with an advanced thermodynamic simulator.The presented method uses this simulator to calculate wellhead pressure, temperature and normalized volume response for a given flow.
An objective of the method used in the AGDs model, is to identify for a given set of input parameters, the "bottlenecks" of the plant which restrict the maximal flow, and to find the device configuration which provides the maximal flow capabilities.The gas flow in or out of the storage may be natural or forced by a set of compressors depending on the pressure difference between the field pipe and the system pipe (Fig. 1).Therefore, the model of the AGDs includes four storage operation modes: injection with compressors, withdrawal with compressors, free flow injection and free flow withdrawal.Depending on the mode, a model of a different set of devices is used to calculate the maximal withdrawal or injection rate.This output value is calculated for a given set of input parameters and device availability which can be imported from the service and maintenance plan.The four main inputs are: field pipe pressure, system pipe pressure, gas mole weight and gas heat capacity.Moreover, in the injection mode case, the gas temperature in the system pipe is needed, while in the withdrawal mode case, the gas temperature in the field pipe is needed.The external input parameters are defined a priori as prediction trajectories, based on some expert knowledge related with an expectation of future gas parameters and system pipe pressure.The length of these trajectories is the same as the time horizon of the simulation.
The field pipe pressure is calculated in each computational step based on information about storage mode, and a set of caverns which are connected to the field pipe.In the withdrawal mode case, the gas temperature in the field pipe is calculated as the weighted average of wellhead temperature of caverns that are connected to the field pipe.
The cavern model is used to calculate the wellhead pressure, temperature and volume response for the given gas flow.Such a model includes the equation of state of gas, gas chemical composition, cavern geometry and the rock mechanical cavern response.The particular, the cavern model used in this simulator was supplied by ESK Gmbh Germany, several other cavern models are also available, such as NOMIX by GreyLogix, CavInfo Software Suite by SOCON, WinKaga by Chemkop, SCGS Toolbox by Technical Toolboxes, Inc.

LIMITATIONS OF CAVERN OPERATION
The main limitations imposed on operation of a single cavern are: -the maximum injection and withdrawal rate, -the range of acceptable cavern pressure, -the maximum allowed pressure change rate during gas injection or withdrawal.These constrains are simple and easy to handle.However, there is also another class of more complex, and possibly time-dependent, cavern operation constraints, which originate from geological restrictions related with rock mechanics, e.g.maximum duration of operation at low pressure levels [4].Violation of such constraints may result in local failures of the rock mass at the cavern wall, leading to the temporary closure of cavern operation, which is an essential loss from the gas storage plant perspective.The detailed form of those restrictions may be different for different plants.One should notice that even a single plant might have caverns associated with different classes of such constraints.The main differences between the classes of caverns are due to their different geometrical properties, such as depth of the cavern location, or the value of v=s ratio.In the example presented in this paper, there is a real life assertion that the plant consists of two classes of caverns.There are 7 class A and 5 class B caverns (the difference between the rock-mechanics restriction for these two cavern classes is discussed in the remainder of this section).The real number of caverns in the plant where the presented method had been tested was different than the number presented here, but the information about the real configuration is protected and cannot be published.The differences and consequences of the rock mechanical recommendations for different cavern properties are discussed in [16].Usage of cavern simulator potentialities in analysis of short-term and long-term mechanical and thermodynamical cavern behavior is discussed in [17].
Figure 2 illustrates real life rock mechanic limitations associated with class A caverns.In the initial stage (before time instance t 1 ) the pressure was in the range p LL ; p max h i .Operating below or at p LL limit causes undesired cavern convergence.Thus reaching or decreasing cavern pressure below this limit is related with a future necessity to maintain the cavern pressure in the range p R ; max h ifor some defined time period.The p LL limit is reached at time t 1 when a dedicated LLUOR cavern timer is started (initial timer state is zero).The system remembers the minimal value of cavern pressure recorded in the time period t 1 ; current time h i .The lowest recorded pressure value determines the maximal allowed time for operating below p LL limit (Fig. 2, orange horizontal lines).In our example, this maximal allowed time is in a range from 100 (when minimal pressure is reached) to 240 days (in case when pressure falls slightly below p LL level).When the LLUOR cavern timer reaches the maximal allowed time, the pressure value must remain in the range of p LL ; p max h i .The nonlinear relationship between the minimum registered pressure and the time limit, is set individually for each cavern and has been approximated by a polynomial of degree sixth.In Figure 2, this time limit, t 0 a ; t 0 b and t 0 c is determined by the cavern state at t a , t b and t c respectively.
For each class A cavern an additional time period (t 2 , t 3 ) is defined.In this time period, the pressure value can be in the range of p LL ; p max h i .At t 2 , the LLUOR cavern timer is frozen.The OFFSET cavern timer starts ticking only if the cavern pressure was outside the restricted range in the last possible day (when LLUOR cavern timer equals the maximal allowed time for staying below p LL level).
It is worth considering the case in which the cavern pressure exits from the pressure range p min ; p LL h ibefore the "maximum allowed time" for staying below p LL has expired.In this case, the timer that keeps track of the operation below p LL is frozen, and all calculated future time points are shifted forward by the length of the time spent above p LL .The OFFSET timer is not ticking, the t 3 time is defined in a such way that the t 3 ; t 4 h i period has some defined length (for example 80 days) and the t 1 ; t 4 h i period is not longer than one year.In the whole time period t 3 ; t 4 h i, the pressure value must be in the range p R ; max h i .All timers are set to zero at time t 4 .The cavern cycle is over and a new cycle is started.When the pressure value reaches the p LL limit for the first time in a new cycle, the LLUOR cavern timer is started, and the whole procedure (that is based on described rules) is repeated.
Figure 3 illustrates real life rock mechanics limitations associated with class B caverns.The cavern pressure is not restricted by geological constraints until it hits the minimal allowed value.At time instance t 1 , the LLUOR timer and the BPR (Bottom Pressure) timer start ticking.The BPR timer counts the days when cavern pressure was at the minimal pressure or close to the minimal pressure.When a cavern reaches minimum pressure the LLUOR timer begins to count the days that cavern pressure is below the p LL level.The maximum value for this timer is defined for each cavern (for example 80 days).The operator can keep the pressure below the p LL level only for a limited time after it hits the minimal pressure.At time t 2 , the cavern pressure becomes higher than minimum pressure plus some margin (), thus the BPR timer is frozen.At t 3 , the pressure exits the restricted range.The LLUOR timer is frozen and the recreation pressure timer starts ticking.From this moment, the pressure has to be in the range p LL ; p max h ifor the time period t 3 ; t 4 h i.The length of this time period is equal to the length of the t 1 ; t 2 h itime period.Thus the recreation pressure timer has to reach the value of the BPR timer.At t 4 , all the timers are set to zero.In the case that the pressure leaves the range min; min þ h iand returns to this range, the length of the t 3 ; t 4 h i time period should be the same as the total time of staying at, or close to the minimal pressure.

CAVERN OPERATION STRATEGY
The goal of this section is to describe a strategy, i.e., the predefined high level decision mechanism for operating caverns.There may be more than one strategy defined for the plant, but only one strategy can be used in a single simulation.For a given state of the gas storage plant, a high level strategy defines a set of actions that can be executed.Each action has a unique priority assigned, which depends on the cost and rock mechanics factors related with the action.In each simulation step, the action which results in the total gas flow that is as close as possible to the gas flow demand is chosen.If there are more desirable actions which result in the same total gas flow, then an action with the highest priority level will be used.Each action has a unique priority level and serves the following three functions: -defines the set of active caverns that can take part in gas injection or withdrawal process; -defines the pressure target for each active cavern; -limits the set of active caverns to those which fit the pressure window p window (Sect.1), the bottom (top) limit of the pressure window is affixed to the lowest (highest) cavern pressure from the set of active caverns in case of injection (withdrawal) mode.
In the optimization algorithm, the set of active caverns (returned by the action) is limited in order to find a minimal set of caverns which provides the highest gas flow.This optimal correction is computed by solving an optimization task (Sect.4) and it is not a part of the strategy definition.In other words, the action of the strategy defines a framework for the optimization task.This framework limits the number of caverns that are used by the optimization task.This limitation takes into consideration plant architecture and physical limits, some expert knowledge related with cavern operation and a mechanism that protects the plant against switching caverns on and off due to caverns thermodynamic effects.

An Example of a Real Life Strategy
In this section, an example of the real life strategy is described.It consists of three substrategies: withdrawal, standard injection and recreation injection substrategy.Each strategy consists of the following two operations: -assigning the discrete system state based on current system parameters, -associating the system state with the set of potential actions by means of the "strategy book", i.e., the rules for decision selection (Tab.1).
In order to handle cavern pressure thermodynamic effects and to avoid oscillations of cavern enable/disable decisions, the following mechanisms were introduced: 1. safety pressure margins on the following parameters: -minimal and maximal cavern pressure, -compressor activation level (please refer to safe compressor activation level in Fig. 4), p LL pressure level (please refer to safe LLUOR level in Fig. 5), p R level, -pressure window size (Sect.1). 2. deadband pressure, which is used to identify an active pressure range in each cavern (Fig. 4 and 5). Figure 6 describes the idea of the deadband pressure.The selected pressure ranges for the withdrawal substrategy for class A and class B caverns (with different rock mechanics rules) are presented in Figure 5 and Figure 4 respectively.The withdrawal substrategy is executed if the gas flow demand for the current step requires the withdrawal mode and no get out action is required.The get out action is required if at least one cavern in the plant requires the injection action in order to fulfill geological constraints related to rock mechanics rules (Sect.2).The set of states and actions has been defined based on selected pressure ranges from Figures 4 and 5.

States for Class A Caverns
State SA 1 .If at least one class A cavern is in the WA 1 pressure range, then the entire set of class A caverns is in the SA 1 state.The following decision is possible in the SA 1 state: -decision DA 1 : select an optimal subset of all class A caverns that are above safe p LL level and get them to safe p LL level.-decision DB 2 : Select an optimal subset of all class B caverns and get them to the safe minimal pressure.Thus the set of active caverns consists entirely of class B caverns.
State identification is based on pressure ranges instead of single pressure levels in order to handle rapid changes of pressure related with thermodynamic effects.Such pressure changes are consequences of gas temperature changes resulting from connecting and disconnecting a cavern from the field pipe.
Table 1 presents the ''strategy book'' which is a set of decision rules.The selection of a decision depends on the state of class A and class B caverns.For instance, if class B caverns are in state SB 2 , and class A caverns are in state SA 1 , then the algorithm will first try to execute the DA 1 action only.If the withdrawal rate for simulated action will be lower than the current demand, then the algorithm will try to use simultaneously both DA 1 and DB 2 actions, i.e., it will select an optimal subset from the set of active caverns that consists of all class A caverns that are above the safe p LL level and all class B caverns.The second combination has a lower priority than using solely the DA 1 decision, but it will be applied if the withdrawal rate provided by this combination is higher.The flow will be spread among selected caverns in proportion to their flow limits.
The standard injection and recreation injection substrategy will be described in a more general form.The framework for defining those substrategies of the default strategy, is the same as for withdrawal mode.The injection substrategy is executed whenever the gas flow demand for a current step requires an injection mode and no get out action is required.The main idea behind this substrategy can be summarized by the following rules: -if there are any class A caverns with pressures in restricted ranges, use such caverns with the highest priority and get them out of restricted ranges; -if there are any class B caverns that reached the minimum pressure and are still below p LL level, use such caverns with the medium priority and get them out of the restricted range; -utilization of all other caverns has low priority.
A combination of the above rules is possible.
The recreation injection substrategy is the last strategy discussed here.It is required if at least one cavern needs to get out from the restricted pressure range in order to fulfill rock mechanics rules.Such a necessity is calculated for each cavern automatically.In the case when the get out action is required, the mode of the plant is set to injection.The plant mode change time (for example 45 minutes) is taken into consideration by the simulation algorithm.This paper does not include a detailed description of this strategy, but the main idea can be summarized by the following rules: High priority rule: if there are any class A or class B caverns that need to get out the restricted range, then set them as active caverns.The pressure target is set in such a way that the algorithm will try to get them out of the restricted range, plus some safety margin that will cover a thermodynamic pressure drop after the get out action is done.This safety margin protects the plant from oscillating between injection and withdrawal modes.
Low priority rule: the set of active caverns consists of all caverns that need to get out of the restricted range and includes all other caverns, starting from the lowest pressure cavern that fits into the allowable pressure window (Sect.1).By construction of the algorithm, the cavern with the lowest pressure from the set of caverns that require a get out action, is always included.However, the cavern gas flow sharing algorithm is different for caverns that require a get out action.The calculation algorithm will try to assign the maximum flow to such caverns, and the remaining flow is spread between all other caverns.The low priority rule was introduced in order to cover the case in which the get out action has been assigned to a small number of caverns, but the injection gas flow demand is high.In such a case, the algorithm might not be able to reach the gas flow demand by applying a high priority rule only.

Simulator Inputs
Main inputs of the simulator are presented below: -scenario of the gas flow demand.One should define a trajectory of withdrawal demand or injection demand for each step of simulation.If the goal is to calculate the maximum possible withdrawal or injection rate on the assumed time horizon (for example 7 days), then this input should be set to infinity in the case of injection mode, or negative infinity in the case of withdrawal mode.To calculate an optimal gas flow for aggregated business nominations, this input should be the same as the aggregated business nominations; -prognosis of physical gas parameters, including system pipe pressure, gas mole weight and gas heat capacity; -scenario of caverns and AGDs availability within the simulation horizon.This information can be taken from the service and maintenance plan; -the initial state of the plant.There is a possibility to set the initial pressure for each cavern manually, or the thermodynamic state of each cavern can be automatically taken from the real plant (there is an assumption that the thermodynamic model of caverns is always synchronized with the real plant).The initial state is also defined by the initial state of rock mechanics rules (Sect.2).One can set the state of each rock mechanics timer manually, or it can be automatically recalculated based on the historical plant data.

Computational Algorithm
The presented method starts the simulation with a particular plant state which has been loaded during the initialization of the algorithm.The plant state comprises the caverns operation history, current configuration and availability of the above ground devices The algorithm chooses the proper action and simulates plant future states with assumed time resolution (e.g. one hour), taking as a further input the gas nomination trajectory.In each step of the main loop the proper action is chosen, the caverns simulation is performed, and the caverns operation history is updated.At the beginning of each step, we identify a discrete plant state according to predefined strategy book.Formally, it is a mapping from the continuous space of variables describing the current pressures in the system, pressure levels of compressor activations and restriction on cavern operations, onto the set of discrete space variables according to which we specify the tentative set of caverns, according to Section 3.1.At this stage, we generate a series of possible action attempts, each attempt gives us a tentative list of ''active'' and ''disabled'' caverns.At this stage, we enter into the subloop where the optimization procedure (Sect.4) is performed for each attempt.At each step of simulation, the model of above ground devices has to be executed several times in order to calculate the maximal possible flow for each probed configuration.The ADG's model contains internal optimization procedures as there are also possible different configurations of devices included in AGD's (for example the serial or parallel connection of compressors).The following work-flow describes in details the procedure in the main loop: A1.Compute the storage mode.First, the algorithm checks if there is a need to start the get out action in order to avoid violating the rock mechanics constraints.If there is such a need then the plant mode is set to "injection".This calculation is based on static pressure-to-volume caverns characteristics and assumed pessimistic estimation of the injection gas flow rate provided by AGDs.If no get out action is required, then the plant mode is determined by the input gas demand trajectory; A2. Compute compressor activation pressure level.
Compressor activation pressure level is the boundary level of gas pressure in the field pipe that is required to start or stop the compressor.It is calculated based on the storage mode and the value of a system pipe gas pressure (storage input pipe); A3.Compute strategy.First, the caverns availability (based on the service and maintenance plan) is checked.Next, the current plant state is matched with the set of states defined by the strategy.By construction, this matching is always possible and unique.The set of possible actions is returned (Tab.1); A4.Simulate each action and choose the one that provides the gas flow that is as close as possible to the gas demand.In case of more than one ''winning'' action, chose the one with the highest priority level.Priorities are set a priori based on economic and safety factors.Let us have a closer look at the A4 procedure.For each simulated action, the following procedures are performed: A4.1.Compute the optimal flow that can be provided by the plant for the given action.This procedure uses the AGDs model and information about cavern flow limits to calculate the maximal possible injection or withdrawal gas flow, and determines an optimal subset of active caverns which provides the highest possible injection or withdrawal gas flow.The optimization task is described in Section 4; A4.2.Compute caverns shares.The share of each selected cavern in the gas injection/withdrawal action is calculated.The set of selected caverns is given by the A4.1 procedure.By default, the shares are proportional to the caverns gas flow limits.The sharing method is different in the case of the get out action.In such a case, our method first wants to use the potential of caverns for which the get out action is set, so it tries to assign the whole gas flow to such caverns.If 100% of the potential of such caverns is used, and still there is some gas flow "unassigned", then the algorithm spreads this flow among all other selected caverns proportionally to their flow limits; A4.3.Assign flows to caverns based on shares and limit the flow demand to the cavern flow limit.If the demand flow assigned to the particular cavern (based on cavern shares) is higher than the flow limit, then it has to be limited.In the default sharing method, the ratio between the flow demand and the flow limit is the same for all caverns.If the limiting of the demand flow occurs, it affects all selected caverns in the same way; A4.4.Simulate the cavern thermodynamic response once the ''winning'' action was chosen.At the initializa-tion of the algorithm, the preceding operational history of all caverns is loaded.At the end of each step, the cavern thermodynamic response of wellhead pressure, volume and temperature for the given flow, is simulated based on the thermodynamic model of the cavern.The cavern states are not only calculated for selected caverns (with assigned flow) but also for all other caverns (zero flow).The cavern operational history is updated after each step of simulation.

OPTIMAL SELECTION OF CAVERNS
The potential maximal flow of a plant is limited by two factors: the sum of gas flow limits of utilized caverns and the possible flow provided by AGDs.The flow provided by AGDs depends mainly on the difference between the system pipe and the field pipe pressure.
The field pipe pressure is dictated by the set of utilized caverns.Thus increasing the number of connected caverns increases the caverns flow potential, but on the other hand increasing it above some boundary value might have a negative effect on the flow provided by the AGDs.In order to determine whether, there is a need to set a "red flag" signaling that the gas should be injected into the cavern initiating the get out action, we assume the worst case scenario, i.e., the minimal gas flow that can be delivered by AGDs and also the possibility that all caverns operating in the restricted pressure range have to simultaneously execute the get out action.The available gas flow is divided among all the caverns performing the get out action, taking into account restrictions such as the maximum allowed pressure change rate during the injection.The general optimization task is defined by finding the following maximum: with constraint: where: x is the set of all selected caverns.This set must at least cover all caverns from an active set that needs to execute a get out action in order to fulfill rock mechanics rules (only if such get out action is required); -f ðxÞ is the maximal possible plant flow as a function of selected caverns; -mðxÞ is the sum of maximal flows provided by selected caverns.The flow provided by each cavern is calculated by choosing a lower value from two possibilities: one is the cavern flow limit, scaled to the size of the simulation step; the second value is the total gas volume which can be still injected or withdrawn from the cavern.For instance, if the gas pressure is close to the minimal allowed pressure in the withdrawal mode, the possible gas flow is low; -p field ðxÞ is the field pipe pressure determined by a set of selected caverns and the plant mode, see Section 1; -g(p field ) is an estimation of the flow limit of AGDs (best case scenario) in function of field pipe pressure.
The AGDs model is used to calculate the maximal possible flow for current p field ; -sðxÞ is the number of caverns in the x set; -c is the injection or withdrawal gas flow demand for the current simulation step; -r is the number of caverns from an active caverns set that needs to execute a get out action in order to fulfill rock mechanics rules; -a is the penalty coefficient; -bðxÞ ¼ 0 if the configuration of utilized caverns is the same as in the previous step, otherwise bðxÞ ¼ 1.
The penalty protects the plant from frequently switching cavern statuses.If the penalty factor is given by a ¼ h a , then the optimizer will change the set of the selected caverns only if such an action will raise the value of the optimized function by more than h a .Otherwise, the optimizer would prefer to stay with the current caverns configuration.By default, h a is equal to the average gas flow rate of a single cavern.
The penalty coefficient can be set to zero according to a penalty reset period which is denoted by s, for example s ¼ 24 hours.Therefore, there is a high probability that the cavern cluster might be changed every s hours.The specific hour is controlled by another variable which is denoted by s offset .
In the first stage of solving, an optimization task of the subset of caverns which can operate in a free flow mode is isolated from the set of all active caverns (the optimizer input).The optimization task is solved separately for this subset of input caverns with compressors turned off, and for the entire set of input caverns with compressors turned on.The action which results in the largest value of the flow is chosen.The value of the function f ðxÞ is calculated in the following way: -compute the maximal flow rate that can be realized by the AGDs for the current mode (withdrawal free flow, withdrawal compressor flow, injection free flow, injection compressor flow), based on: -the calculated decision about compressor utilization, -the information about AGDs availability (i.e. from the service and maintenance plan), -the assumed trajectory of gas physical parameters (i.e.gas pressure in a transport line), -the computed field pipe temperature (based on wellhead gas temperature from the set of active caverns) in case of withdrawal mode; -scale the computed flow rate to the simulation time resolution; -if the maximal possible flow is higher than the gas flow demand for the current simulation step, then the maximal gas flow is reduced to the assumed gas demand.If there is a get out action required then the maximal gas flow provided by the plant is calculated in the following way: -let u be a sum of all injection flow limits related with caverns that need to get out from restricted ranges (according to rock mechanics limitations); -let d be a gas injection demand for a current step; -let z be a maximal gas flow provided by the AGDs; -final flow is a min(z,max(u,d)).
The calculation of an optimal subset of active caverns can take a long time.To test each subset of active caverns, the algorithm has to execute the model of the AGDs.However, in many cases we can assume that if the difference between field pipe pressure and system pipe rises, the withdrawal flow rises and the injection flow drops.With such an assumption the optimization task can be solved in the following way: -calculate the flow of the plant for the set of selected caverns in the previous simulation step (only if such set is allowed by strategy and rock mechanics constraints).In such a case, there is no penalty; -sort all the active caverns in such a way that the first element of the sorted list is a cavern with the lowest (highest) pressure in case of an injection (withdrawal) mode; -in the first iteration calculate, the function f ðxÞ, where x consists of the first element of the sorted list; -in all subsequent iterations of the optimization task: -add the next cavern from the sorted list to the set x, -if the set x is the same as the set of selected caverns from the previous step, skip this iteration ("no penalty" solution was already computed); -iterations are stopped if the value of the function f ðxÞ does not increase in the next iteration.If the increase is lower than the adjustable parameter d the iterations are also stopped.Please note that if the penalty value is constant (the set of selected caverns form the previous simulation step is excluded from calculations), the function f ðxÞ is convex under the presented assumptions.Thus the calculation of an optimal subset of caverns can be even speeded up in case of a large caverns set by applying one of many discrete optimization methods dedicated to convex functions.
The presented algorithm calculates the optimal subset of active caverns that provides the highest flow.In case of more than one ''winning'' subset, the algorithm will choose the subset with the lowest number of caverns.

SIMULATION RESULTS
The long term (three months) simulation of the plant model was executed in order to show the method decisions and trajectories of the plant parameters.The step size was 1 hour.In the presented illustrative example, gas flow nominations were set in the following way: -in steps 0-39, the withdrawal flow demand was 2.5 million normal m 3 =h, which exceeded the maximal plant flow, therefore, the actual flow was determined by the optimization procedure; -in steps 40-838, the withdrawal flow demand was between 0.2-0.5 million normal m 3 =h.Supplying this demand was within the capabilities of the plant; -in steps 839-1 456, gas injection was required, the nominations were 0.2-0.5 million normal m 3 =h.Such an injection was within the capabilities of the plant; -in steps 1 457-2 160, gas injection of 2.5 million normal m 3 =h was required, which exceeded the maximal flow of the plant, therefore the actual flow was determined by the optimization procedure.
For the purpose of demonstrative interpretation of simulator actions, the gas pressure in the system line was kept constant.In the initial state, the amount of gas in each cavern was sufficiently large to allow the cavern pressure to be well above the p LL level (Fig. 2 and Fig. 3).The resulting pressure trajectories are shown in Figure 7.
In the first stage, the gas free flow is possible (withdrawal without the need to turn compressors on).The optimizer procedure chooses the set of caverns that provides the maximal flow of the plant.Note that not all available caverns are used, the flow is restricted by the AGDs operational limits (Sect.5).When the differences between the pressures in the cavern pool become large enough to induce the switch of the configuration, the caverns with the lowest pressure are turned off and the caverns that were inactive are turned on.Around step 100, the cavern pressure values approach the compressor activation level.The optimizer performs more frequent switches in order to maintain the configuration that allows the free flow.After step 229 (point a in Fig. 7), free flow is no longer possible and the compressor has to be turned on.From step 259 (point b in Fig. 7), only class B caverns are used, which is dictated by high priority action for states SA 2 and SB 2 according to the ''strategy book'' (Tab.1).After step 321 (point c in Fig. 7), the current configuration does not provide the demanded flow, therefore the set of class A caverns is selected (medium priority action) which provides the demanded flow.At step 421 (point d in Fig. 7), the flow demand has decreased from 0.45 to 0.2 million normal m 3 =h, therefore the high  The magnification of pressure trajectories during the injection action.The configuration switching occurs every s (24 hours) which is denoted by stars on x-axis.
priority action is chosen again as it covers 100% of the flow demand.This action is kept until step 543 (point e in Fig. 7), when the medium priority action is chosen again.At step 771 (point f in Fig. 7), one of the class B caverns starts the get out action, as is dictated by the state of the LLUOR timer and the current pressure in the cavern.The plant mode is then changed to injection in spite of the withdrawal demand.The entire flow is directed into the cavern performing the get out action.From step 839 (point g in Fig. 7), gas injection is required (the gas demand varies from 0.2 to 0.5 million normal m 3 =h), therefore additional caverns are connected.At step 882 (point h in Fig. 7), the get out action is completed and the plant continues to fulfill the injection requirements.At step 1 457 (point i in Fig. 7), the injection demand rises above the plant's available potential, therefore the optimal configuration to achieve the maximal plant flow is used.The injection into each cavern takes place until the pressure in the cavern reaches the maximal allowable pressure value, with configuration switching occurring every s hours, where s ¼ 24 (Fig. 8). Figure 9 illustrates the improvement of the plant's flow when the optimal configuration of caverns is used instead of all available caverns connected.This improvement depends obviously on several different parameters, we present it as a function of the pressure spread in the caverns set (represented by pressure standard deviation).The mean value of the pressure values was kept constant.The flow's improvement due to the optimizer action rises as a function of the pressure spread in the set of active caverns up to when the standard deviation of the spread was about 15 bar and the pressure values cover the entire allowable pressure window.

CONCLUSIONS
The presented method constitutes a framework which bridges models of plant elements with the operational rules.The main achievement is an integration of different computational layers into an independent fully functional tool.The presented method defines an optimal strategy (decision rules) for a plant operation.Such an optimal strategy significantly increases the gas storage potential without any change to infrastructure.The increased capacity results in more dynamical responses to changes in natural gas demand, and hence the gas availability for the end user.This method allows for computation of realistic injection and withdrawal plant flow limits that can be transformed to the storage user's nomination limits.Moreover, SOE can be used for calculating an optimal realization of aggregated business nominations for the plant.Proper simulation and identification of the plant's potential is only possible if we consider important nonlinear limits related with cavern operation and the optimal strategy (set of decision rules) which controls the plant operation.It should be pointed out that, in comparison to other gas storage plant simulators, the SOE method takes into account time-dependent recommendations, resulting from rock mechanics, that are related to gas cavern operation.These rules vary for different types of caverns.SOE monitors the state of each cavern following these recommendations.The simulator automatically detects the necessity to execute an injection that will protect a given cavern from violating geological and operational constraints, therefore SOE takes responsibility for meeting the complex nonlinear geological constraints.Optimizing the operation of depleted or partially depleted gas field reservoirs is in principle a different issue [18], however the general framework presented here might be adapted also for that case.The SOE method has been initially tested in an underground gas storage plant located in Germany where it assists the daily operation.The gain in the maximum flow (between optimized and non-optimized configuration) depends significantly on the pressure spread in the caverns set.
Figure 1Simplified model of the plant.
Figure 2Rock mechanics limitations associated with class A caverns.

Figure 3
Figure 3 Rock mechanics limitations associated with class B caverns.

Figure 4 Figure 5
Figure 4 Default strategy, withdrawal: selected pressure ranges for B class caverns.

Figure 7
Figure 7Pressure trajectories of class A caverns (green solid line) and class B caverns (red dashed line).The explanation of performed actions (a-h) and specifications of different gas demand nominations (rectangles above x-axis and vertical orange lines) are presented in Section 5.

TABLE 1
The "strategy book".Decision (or decision combination) selection depends on states of class A and class B caverns Thus the entire set of active caverns consists of class A caverns that are above safe p LL level.State SA 2 .If all class A caverns are in the WA 2 or WA 3 range, then the set of class A caverns is in the SA 2 state.The following decision is possible in the SA 2 state: -decision DA 2 : select an optimal subset of all class A caverns and get them to safe minimal pressure.Thus the set of active caverns consists entirely of class A caverns.State SB 1 .If at least one class B cavern is in the pressure range, then the entire set of class B caverns is in the state.The following decisions are possible in the SB 1 state: -decision DB 1 : select an optimal subset of all class B caverns that are above safe compressor activation level and get them to the safe compressor activation level.Thus the set of active caverns consists of class B caverns that are above safe compressor activation level; -decision DB 2 : select an optimal subset of all class B caverns and get them to the safe minimal pressure.Thus the set of active caverns consists entirely of class B caverns.State SB 2 .If all class B caverns are in the WB 2 or WB 3 range, then the set of class B caverns is in the SB 2 state.The following decision is possible in the SB 2 state: